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ABSTRACT 

The  Distributed  Array  Processor  (DAP)  manufactured  by  International 
Computers  Limited  is  an  array  of  1-bit  2 00- nanosecond  processors*  The  Pilot 
DAP  on  which  the  present  work  was  done  is  a  32  x  32  array;  the  coemerclally 
available  machine  is  a  64  x  64  array.  He  show  how  the  projected  SOR  algorithm 
for  the  linear  complementarity  problem  Aw  b,  w  £  0,  wT(Aw  -  b)  ■  0,  can  be 
adapted  for  use  on  the  DAP  when  A  is  the  'finite-difference'  matrix 
corresponding  to  the  difference  approximation  to  the  Laplace  operator. 
Application  is  made  to  two  linear  complementarity  problems  arising,  respec¬ 
tively,  from  two- and  three-dimensional  porous  flow  free  boundary  problesis. 
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SIGNIFICANCE  AND  EXPLANATION 


An  array  processor,  the  Distributed  Array  Processor  manufactured  by 
International  Computers  Limited,  has  recently  become  available#  The  DAP  is  an 
array  of  1-bit  processors;  in  the  production  machine  there  are  4096  processors 
arranged  as  a  64  x  64  array#  It  is  normally  programmed  in  an  array  processing 
extension  of  Fortran. 


It  is  of  interest  to  develop  algorithms  which  can  efficiently  use  the 
great  confuting  capacity  of  the  DAP.  We  show  how  the  projected  SOR  algorithm 
for  the  linear  complementarity  problem  Aw  >_  b,  w  >_  0,  wT(Aw  -  b)  ■  0  can  be 
adapted  for  use  on  the  DAP  when  A  is  the  finite- difference  matrix 
corresponding  to  the  difference  approximation  to  the  Laplace  operator. 
Application  is  made  to  two  linear  complementarity  problems  arising, 
respectively,  from  two  and  three-dimensional  porous  flow  free  boundary 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  me,  and  not  with  the  authors  of  this  report. 


THE  SOLUTION  OF  LINEAR  COMPLEMENTARITY  PROBLEMS  ON  AN  ARRAY  PROCESSOR 


C.  W.  Cryer* ' 1 ' 2,  P.  M*  Flanders'1-,  D.  J.  Hunt+,  S.  F.  Reddaway*, 

and  J.  Stansbury**' ^ 

1*  Introduction 

An  LCP  ( linear  complementarity  problem)  is  a  problem  o£  the  form:  Find  a  real 
n-vector  w  »  (wi)  satisfying 

(a)  Aw  _>  b  , 

(b)  w  >  0  ,  (1.1) 

(c)  wT(Aw  -  b)  -  0  , 

where  b  -  (b^)  is  a  known  real  n-vector  and  A  “  (a^j)  is  a  known  real  n  x  n 
matrix. 

Linear  complementarity  problems  arise  in  many  contexts  (Balinski  and  Cottle 
[1978]).  In  particular,  there  is  a  close  connection  between  linear  complementarity 
problems  and  variational  inequalities  (Cottle,  Giannessi  and  Lions  [1980],  Cryer  and 
Dempster  [  1 980 ] ) . 

Many  problems  in  continuum  mechanics  can  be  reformulated  as  variational 
inequalities  (Duvaut  and  Lions  [1976],  Kinderlehrer  and  Stampacchia  [1980]),  which, 
when  discretized,  reduce  to  linear  complementarity  problems  of  the  form  (1.1)  with 
special  features! 
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1*  A  is  a  large  matrix,  perhaps  of  order  25,000. 

2.  A  is  a  'finite-difference'  or  'finite-element'  matrix;  in  particular, 

A  is  sparse  with  a  great  deal  of  structure.  (1*2) 

3.  A  large  percentage  of  the  elements  of  the  solution  w  are  non-zero. 

Because  of  these  special  features,  the  standard  methods  of  solving  linear 
complementarity  problems  are  not  very  efficient,  and  methods  of  solution  have  been 
developed  which  take  advantage  of  the  structure  of  A:  projected  SOR  (Cryer  [1971], 
Glowinski  [1971]);  modified  block  SOR  (Cottle,  Golub,  and  Sacher  [1978]);  multigrid 
projection  (Brandt  and  Cryer  [1980]);  and  generalizations  of  projected  SOR 
(Mangasarian  [1977]).  Cryer  [1980]  briefly  surveys  much  of  this  work. 

In  the  present  paper  we  consider  the  use  of  the  parallel  computer  DAP  to  solve 
linear  complementarity  problems  with  the  features  (1.2).  The  DAP  (Distributed  Array 
Processor,  manufactured  by  International  Computers  Limited) ,  which  is  an  SIMD  array 
of  typically  64  x  64  processors,  is  described  in  Section  2.  In  Section  3  we 
describe  the  implementation  on  the  DAP  of  projected  SOR  to  solve  a  linear 
complementarity  problem  derived  from  a  two-dimensional  porous  flow  free  boundary 
problem,  and  in  section  4  we  extend  this  work  and  solve  a  linear  complementarity 
problem  derived  from  a  three-dimensional  porous  flow  free  boundary  problem.  In 
Section  5  we  comment  on  possible  future  developments,  and  the  overall  conclusions 


are  in  Section  6 
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2*  The  Pilot  DAP  (Distributed  Array  Processor) 

The  preeent  work  was  carried  out  on  the  Pilot  32  x  32  DAP  at  Stevenage# 

England,  and  we  will  describe  this  machine  first.  A  64  x  64  version  is  available, 
and  the  minor  differences  between  the  two  machines  are  indicated  at  the  end  of  this 
section. 

PAP  Hardware 

The  essential  features  of  the  Pilot  DAP  hardware  are  as  follows  (Flanders, 

Hunt,  Reddaway,  and  Parkinson  11977],  Reddaway  [1979]): 

1.  A  32  x  32  array  of  identical  processing  elements  (PEs)  with  a  cycle  time  of  200 
nanoseconds . 

2.  Each  PE  has  a  one-bit  adder,  2K  bits  of  storage,  and  three  one-bit  registers  (a 
general  purpose  register  for  accessing  data  and  performing  arithmetic;  a  carry 
register;  and  an  activity  control  register). 

3.  Bach  PE  is  connected  to  its  four  neighboring  PE's  (North,  South,  East,  and 
West).  In  a  given  cycle  all  PE's  access  their  neighbor  in  the  same  direction 
(determined  by  the  program).  In  addition,  the  PEs  are  linked  by  row  and  column 
highways  which  connect  together  all  the  PEs  in  each  row  and  column. 

4.  There  is  a  master  control  unit  (M 03)  which  broadcasts  instructions  to  all  the 
PEs.  All  PEs  can  perform  the  same  instruction  simultaneously,  but  certain 
instructions  are  only  effective  if  the  activity  control  register  is  ' true' . 

PAP  Software 

A  programm  to  run  on  a  PAP  system  normally  comprises  a  standard  FORTRAN  program 
and  a  number  of  subroutines  and  functions  written  in  an  array  processing  extension 
at  FORTRAN  known  as  DAP-FORTRAN  (Flanders  [1979],  Gostick  [1979],  ICL  [1979]). 

The  standard  FORTRAN  is  executed  by  the  host  computer  and  provides  mainly 
input-output  and  overall  control.  The  PAP-FORTRAN  is  executed  by  the  DAP  and 
provides  high  speed  computation.  Data  is  shared  between  them  using  common  blocks 


held  in  DAP  store 


Some  features  of  DAP-FORTRAN  are  described  below 


In  addition  to  the  data  types  of  FORTRAN ,  DAP-FORTRAN  has  two  new  data  types t 
vector  and  Matrix.  With  a  32  *  32  DAP  a  vector  has  32  components  and  a  matrix 
has  32  x  32  components*  the  components  can  be  real,  integer,  or  logical. 

For  example,  the  data  statements 

REAL  U(  ),  V(  ,  ),  W( ,5) ,  X( , ,3 ) , 

INTEGER  A( , 1 ) ,  B(  ),  C(,,4)  (2.1) 

LOGICAL  FLAGS ( , 2 ) ,  MASK(  ,  ) 

declare  D  (a  real  vector),  V  (a  real  matrix),  W  (an  array  of  5  real 
vectors),  X  (an  array  of  3  real  matrices),  A  (an  array  of  1  integer  vector),  B 
(an  integer  vector),  C  (an  array  of  4  .integer  matrices),  FLAGS  (an  array  of  2 
logical  vectors),  and  MASX  (a  logical  matrix). 

Expressions  in  DAP-FORTRAN  can  consist  of  scalars,  vectors,  and  matrices  with 
the  usual  unary  and  binary  operations.  Operations  on  vectors  and  matrices  are 
performed  in  parallel  using  all  32  x  32  pbs. 

Operations  between  a  scalar  and  a  vector  or  a  matrix  cause  implicit  expansion 
of  the  scalar  to  the  necessary  dimensions.  For  example,  if  M  is  a  matrix  of  size 
32  x  32  and  S  is  a  scalar,  then  M  ■  M  +  S  causes  S  to  be  implicitly  expanded 
to  size  32  x  32  with  each  element  being  equal  to  S;  then  the  corresponding  elements 
of  "matrix"  S  and  matrix  M  are  added  in  parallel  and  assigned  to  M  in  parallel. 

Arrays  of  vectors  and  matrices  may  be  used  to  construct  more  complex 
structures.  To  process  a  vector  or  matrix  array  requires  performing  calculations  on 
the  individual  vectors  or  matrices  in  the  array. 

Selection  and  updating  of  parts  of  vectors  and  matrices  can  be  performed  using 
the  powerful  indexing  capabilities  of  DAP-FORTRAN.  Matrix  sections  can  be  specified 
by  emitting  subscripts  along  which  all  elements  are  to  be  taken.  Using  this,  whole 
rows  or  coltams  can  be  selected  from  matrices.  For  example,  M(I,)  specifies  the 


I-th  row  of  matrix  M 


Shift  indexing  is  a  very  useful  feature  of  DAP-FORTRAN*  For  example,  in  a 
simple  solution  of  Laplace's  equation  on  s  32  v  32  grid  we  wish  to  replace  each 
element  with  the  average  of  its  four  neighbors*  This  could  be  coded  in  FORTRAN  as: 


DO 

10 

I 

-  2,31 

DO 

10 

J 

“  2,31 

Y(I,J)  -  (X(I  +  1,J)  +  X(I  -  1,J)  +  X(I,J  +  1)  +  X(I,J  -  1))  /  4.0 
10  continue 

Further  code  would  be  needed  to  handle  elements  on  the  edges  of  the  matrix* 

The  DAP-FORTRAN  code  is  much  simpler: 

X  -  (X(+,)  +  X(-,)  +  X(,+)  +  X(,-))  /  4.0  .  (2.1) 

The  term  X(+,)  uses  shift  indexing.  In  particular,  X(+,)  specifies  a  matrix 
where  the  11,  J)  element  is  the  (I  +  1,J)  element  of  X,  for  1  £  I  £  32  and 
1  <  J  £  32.  Thus,  X(+,)  contains  all  the  "south"  neighbors  of  X.  Edge  values 
( corresponding  to  subscripts  0  or  33)  are  defined  to  be  zero.  As  an 
alternative,  cyclic  geometry  may  be  specified  by  using  a  GEOMETRY  statement. 

Longer  shifts  can  be  performed  by  explicit  system  functions;  for  example, 
SHS(X,I)  shifts  the  matrix  X  I  positions  to  the  south.  Note  that  since  all  the 
updating  is  performed  simultaneously,  it  is  not  necessary  to  write  the  results  to 
another  matrix. 

Logical  matrices  and  vectors  can  be  used  to  select  elements  from  an  array.  For 
example,  if  we  wished  to  update  only  certain  elements  of  X  in  statement  (2.1),  we 
could  set  the  corresponding  elements  of  LM,  a  logical  matrix,  to  true  and  all 
other  elements  of  LM  to  false.  That  is,  if  X(I,J)  is  to  contain  the  average  of 
its  four  neighbors,  then  LM(I,J)  is  set  to  true.  Otherwise,  LM(I,J)  is  false. 
Then  the  following  statement  performs  the  required  task: 

X(LM)  -  (X(+,)  +  X(-,)  +  X(,+)  +  X(,-))  /  4.0  . 
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DAP-FORTRAN  has  a  number  of  useful  system  functions  whose  arguments  and  results 
may  be  scalars,  vectors,  or  matrices.  The  ALTC,  ALTR,  MERGE,  MAX,  and  ABS  func¬ 
tions  will  be  briefly  described  since  these  are  used  in  the  programs  in  this  paper. 

The  functions  ALTC  and  ALTR  return  logical  matrices.  If  C  is  the 
argument  to  ALTC,  then  the  first  C  columns  of  the  result  matrix  are  set  to 
false,  the  next  C  columns  to  true,  the  next  C  columns  to  false,  etc.  ALTR 
performs  similarly  for  rows. 

The  function  MAX  (now  named  MAXV)  returns  a  scalar  equal  to  the  largest 
number  in  its  vector  or  matrix  argument.  The  function  ABS  returns  a  vector  or 
matrix  containing  the  absolute  value  of  every  element  in  its  argument. 

The  function  MERGE  takes  three  arguments  and  returns  a  matrix.  The  first  two 
arguments  are  matrices  (or  implicitly  expanded  scalars)  and  the  third  argument  is  a 
logical  matrix.  If  the  (I,J)  element  of  the  logical  matrix  is  true  then  the 
(I,J)  element  of  the  result  matrix  is  set  equal  to  the  (I,J)  element  of  the  first 
matrix;  otherwise,  it  is  set  equal  to  the  (I,J)  element  of  the  second  matrix. 

Examples  of  DAP-FORTRAN  programs  are  given  in  Sections  3  and  4,  and  the 
Appendices . 

DAP  Arithmetic 

When  a  DAP-FORTRAN  program  is  executed  by  the  DAP,  expressions  involving  only 
scalars  are  executed  sequentially,  but  operations  on  vectors  and  matrices  are 
performed  in  parallel  by  the  PEs. 

The  DAP  memory  can  be  visualized  as  a  cuboid,  with  2K  horizontal  planes,  each 
plane  being  a  32  x  32  square  of  bits.  The  32  x  32  array  of  PEs  lies  on  top  of  the 
cube,  and  each  column  of  2K  bits  belongs  to  the  PE  above  it. 

Two  storage  modes  are  used  in  DAP-FORTRAN:  vertical  and  horizontal.  Scalars 
and  vectors  are  stored  in  horizontal  mode  while  matrices  are  held  in  vertical  mode. 

In  vertical  mode,  each  number  is  held  entirely  within  the  store  of  one  PE  with 
successive  bits  in  successive  store  locations.  Thus,  for  an  integer  matrix,  the 


sign  bit  of  every  element  in  the  matrix  would  be  held  in  the  same  store  address  of 
each  PE. 

Xn  horizontal  mode,  a  number  is  spread  along  a  row  of  PEs.  Thus,  a  scalar 
occupies  one  row  while  a  vector  occupies  32  rows.  DAP  instructions  are  also  stored 
in  this  format. 

All  arithmetic  is  carried  out  using  subroutines.  Some  operation  times  for  32 
bit  numbers  are  given  in  Table  2.1. 

It  will  be  noted  that  vector  arithmetic  is  faster  than  matrix  arithmetic.  This 
is  because  a  row  of  PEs  are  available  for  each  vector  component,  while  only  one  PE 
is  available  for  each  matrix  component. 

Some  of  the  quoted  computation  times  are  data  dependent.  In  particular,  matrix 
multiplication  by  a  scalar  typically  varies  from  170)s  to  200)18  depending  upon 
the  distribution  of  zeros  in  the  binary  representation  of  the  constant;  for  special 
scalars  such  as  .5  or  3  the  multiplication  time  can  be  as  low  as  60)s. 


Operation 

Matrix 

Vector 

Scalar 

floating  point  addition 

140-180)18 

54)  s 

27)s 

floating  point  multiplication 

315)s 

50)s 

34)s 

floating  point  multiplication 
by  a  scalar 

60-200us 

40)  s 

- 

One  shift  of  a  real  matrix, 
e.g.  X  ( + , ) « 

15|is 

2)8 

- 

Move  a  floating  point  matrix 

15)s 

2)s 

2)s 

logical  AND 

2)8 

2)s 

2)s 

logical  mask 

1)S 

2)s 

- 

Note:  Times  are  slightly  different  on  production  DAPs. 

Table  2.1.  Average  DAP-FORTRAN  arithmetic  times  for  the  Pilot  DAP. 
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Host-DAP  Interface 

The  sequence  o£  operations  for  compiling  and  running  DAP  programs  is  as 
follows. 

(a)  The  host  computer  compiles  the  host  FORTRAN  program  and  the  DAP-FORTRAN 
subroutines  into  host  and  DAP  machine  codes  respectively. 

(b)  DAP  machine  code,  incorporating  all  necessary  low  level  subroutines,  is  loaded 
into  DAP  memory  in  horizontal  mode  where  it  occupies  a  few  bits  of  each  PE's 
memory.  Host  machine  code  is  loaded  into  the  host  memory. 

( c)  Execution  begins  in  the  host  and  control  is  transferred  to  the  DAP  as  required 
by  subroutine  calls.  On  completion  of  DAP  processing  the  host  resumes 
execution  at  the  point  following  the  call. 

Detailed  information  on  the  Pilot  DAP  relevant  to  understanding  the  programs  in 
this  paper  is  given  in  Appendix  A. 

The  Production  DAP 

The  current  production  DAP  is  generally  similar  to  the  Pilot  but  differs  as 
follows : 

(a)  there  are  4096  PEs  arranged  in  a  64  x  64  array? 

(b)  each  PE  has  4K  bits  of  memory? 

(c)  arithmetic  operations  differ  somewhat  in  timing  but  are  overall  a  little 
faster? 

(d)  coupling  between  host  and  DAP  is  more  direct  so  the  interface  is  simpler  than 
indicated  in  Appendix  A. 
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3.  Numerical  solution  of  a  tw^-dimensional  free  boundary  problem 

The  flow  of  water  through  a  porous  dam  is  a  well-known  model  problem, 
seeps  from  a  reservoir  of  height  H  through  a  rectangular  dam  of  width  L 
reservoir  of  height  h.  Part  of  the  dam  is  saturated  and  the  remainder  of 
is  dry.  The  wet  and  dry  regions  are  separated  by  an  unknown  free  boundary 
must  be  found  as  part  of  the  solution. 


As(0,H) 


FslL.H) 


-  HEAD  — 

-  WATER -- 


^.1 


FIG.  3.1.  Flow  through  a  porous  rectangular  dam  R. 


i 
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As  shown  by  Baiocchi  [1972]  the  problem  can  be  formulated  as  follows: 
u  on  the  rectangle  R  =  ABCF  such  that 

(a)  -V2u  _>  -1 ,  on  R  , 

(b)  u  >_  0,  on  R  , 

(c)  u(-V2u  +  1)  =  0,  on  R  , 

and 

(H  -  y)2/2  on  AB  , 

(h  -  y)2/2  on  CD  , 


[H2(L  -  x)  +  h2x] /2L,  on  BC  , 
0,  on  DFA  . 


Water 
to  a 
the  dam 
T  which 


Find 

(3.1) 


(3.2) 
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The  wet  region  of  the  dam  consists  of  the  points  where  u  >  0  and  the  dry  region 
consists  of  the  points  where  u  ■  0. 


When  the  problem  (3.1),  (3.2)  is  approximated  using  the  classical  five  point 
difference  approximation  for  the  Laplace  operator,  one  obtains  an  LCP  of  the  form 
(1.1),  where  the  matrix  A  and  right  hand  side  b  are  the  same  as  those  that  would 
be  obtained  if  the  Dirichlet  problem 


-V  u  =  -1 ,  on  R  , 


u  =  g,  on  R  , 


(3.3) 


(3.4).  One  generates  a  sequence  of  approximations  w 


(a) 


w 


<k) 

i“1*j 


+  w 


<k) 

i+1/j 


+  w 


(k) 

i/j-1 


+  w 


<k) 

i,j+1 


(Ax)2  , 


(b) 


w<)c+1^) 

ij 
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(c)  wij+1 }  “  max(°»w[j+  ^  >>  *  (3.5) 

for  1  <  i  <  M  and  1  <  j  <  N  » 

(d)  wij*1}  "  9i;j/  for  ((j  -  1  )Ax, ( i  -  1)Ax)  e  3r  . 

It  is  known  that  the  projected  Jacobi  method  will  converge  (Mangasarian  [1977]). 

If  M  <  32  and  N  <  32  then  the  gridpoints  can  be  regarded  as  a  subset  of 
a  32  x  32  array,  and  one  PE  can  be  associated  with  each  gridpoint.  Defining  w*k*, 
w(k+1)  and  z(k)  a8  reai  DAP-FORTRAN  matrices,  the  computation  (3.5)  is  trivial  to 
implement  on  the  DAP. 

In  Figure  3.2  we  list  a  DAP  subroutine  JACOBI  which  solves  the  dam  problem  for 
the  case  h«  0,  H~31,  L*31,  M«N«  32,  and  Ax  -  1.  This  subroutine  could  be 
called  by  a  host  program,  which  could  then  print  the  answers  in  the  matrix  W. 

Using  the  operation  times  given  in  Table  2.1  we  can  readily  estimate  the  time 
required  per  iteration  in  the  main  loop  of  the  JACOBI  subroutine  (see  Figure  3.3). 
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SUBROUTINE 


JACOBI 


LOGICAL  MASK(  ,  ),  WSIGN(  ,  ) 
REAL  W(  ,  )  ,  Z(  ,  ) 

REAL  INDEX (  ) 

EQUIVALENCE  (W,WSIGN) 

HEIGHT  -  31.0 

WIDTH  *  31.0 

DO  10  I  -  1,32 

INDEX  ( I )  »  (32  -  D/31.0 

10  CONTINUE 

W  -  0 

TEMP  *  HEIGHT*HEIGHT* . 5 
W( 1 ,  )  -  TEMP  *  INDEX 

W(  ,1)  «  TEMP  *  INDEX  *  INDEX 
MASK  =  .TRUE. 

MASK( 1 ,  )  *  .FALSE. 

MASK (32,  )  =  .FALSE. 

MASK(  , 1)  -  .FALSE. 

MASK(  ,32)  -  .FALSE. 

DO  50  1*1,  100 

1  Z  *  W(+,  )  +  W(-,)  +  W(  ,+) 

+  W(  ,-)  -  1.0 

2  W( MASK)  *  .25*Z 

3  W(MASK  .AND.  WSIGN)  -  0.0 
50  CONTINUE 

END 


Declare  logical  32  x  32  matrices,  MASK  and 
WSIGN 

Declare  real  floating  point  32  x  32  matrices 
W  and  Z 

Declare  a  real  floating  point  32-vector  INDEX. 
Declare  the  logical  matrix  WSIGN  equivalent  to 
the  first  bit,  the  sign  bit,  of  the  matrix  W. 


Initialize  INDEX  vector. 


Clear  matrix  W 


Set  values  of  the  matrix  W  equal  to  g 
on  bottom  ( BC ) . 

Set  values  of  the  matrix  w  equal  to  g 
on  left  (AB) . 


Set  the  matrix  MASK  to  be  true  at  interior 
points  and  false  at  boundary  points. 


Start  of  main  loop 

Sum  neighbors  and  store  in  Z  matrix. 

Transfer  average  to  W  at  interior  points. 

Project  by  setting  W  *  0  at  points  where 
MASK  is  true  and  the  sign  of  W  is  negative. 


FIG.  3.2.  The  DAP  subroutine  JACOBI 


Statement 


Operations 


Time  (us) 


Z  -  W(+,  )  +  W(-,  )  +  W(  ,+) 

+  W(  ,-)  -  1.0 


4  floating  point  matrix 
additions/ subtractions 
4  index  shifts 
1  scalar-matrix  assignment 


640 

60 

15 


W(MASK)  -  ,25*Z 


1  floating  point  matrix  multipli¬ 
cation  by  a  special  constant 
1  logical  mask 


70 

1 


W(MASK  .AMD.  WSIGM )  »  0.0 


1  logical  AMD 
1  logical  mask 
1  scalar-matrix  assignment 


2 

1 

15 


DO  50  I  -  1,100 


_ 7 

811 


PIG.  3.3.  Estimated  computation  time  for  the  main  loop,  of  JACOBI. 


From  Figure  3.3  we  see  that  one  projected  Jacobi  iteration  over  the  whole 
32  x  32  grid  requires  Slips. 


The  projected  SOR  method 


Let  w*°*  «  (wj°* )  be  an  initial  guess  for  the  solution  w  *  (wij)  of 


ij 


(3.4).  In  the  usual  implementation  of  projected  SOR  one  generates  a  sequence  of 


.  .  .  (k)  (k)  _  ,, 

approximations  w  *  wij  as  follows t 


(a) 


(k)  (k+1 )  (k)  ^(k+1)  .  (k)  ..  .2 

ij  i-1,j  i+1,j  i,j-1  i,j+1  “  (Ax)  ' 


(b) 


C^V2)  .  w(k)  +uj(z(k)  .  4w(k)}  /4  f 


ij 


ij 


-  (w^Jz**1  +  (1  -  uW***  ' 


(3 


(c) 


m  maxtO.w^  ^ 


for  1  <  i  <  M  and  1  <  j  <  N  , 


where  u  is  a  constant,  the  over-relaxation  parameter 


It  is  known  that  the  iteration  (3*6)  converges  for  all  initial  guesses  w^°  ^ 

iff  0  <  u  <  2  (Cryer  [1971],  Glowinskl  [1971]). 

The  implementation  (3.6)  is  not  suitable  for  parallel  computation  because  the 

new  values  w*k+1  ^  cannot  be  computed  simultaneously!  and  wj***1]  must  be 

i“i,3  i,3“i 

known  before  wj*+1 '  can  be  computed. 

However,  there  is  a  simple  but  ingenious  way  of  making  SOR  suitable  for 
parallel  computation.  In  the  implementation  (3.6),  we  order  the  gridpoints  by  rows 
and  columns  (Figure  3.4a).  Instead,  let  us  visualise  the  gridpoints  as  forming  a 
red-black  chess  board  and  number  first  the  red  points  and  then  the  black  points 
(Figure  3.4b). 


13  14  15  16 


9  10  11  12 


5  6/8 


12  3  4 


fisl  G)  Rel  (a) 

©  DU  ©  Q±) 

B3  ®  GH  © 

©  0  ®  mo 


(a)  Usual  (b)  Red  and  Black  1  1 

FIG.  3.4.  Orderings  of  gridpoints  ( for  a  4  x  4  grid) . 


Applying  projected  SOR  to  the  points  numbered  as  in  Figure  3.4(b)  we  find  that 
each  projected  SOR  Iteration  can  be  broken  down  into  two  stages!  in  the  red( first) 
stage  projected  SOR  is  applied  to  the  red  points!  and  in  the  black( second)  stage 
projected  SOR  is  applied  to  the  black  points i 
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Rad  8ta< 


(k,red)  (k# black)  (k,black)  (k,blaok)  (k,black) 

ij  i+1,j  1-1. j  i,j+1  i.j-1 


(4k)2  , 


(b)  'r*d'  -  (»/4)«;|;'t•'4,  *  (i  -  »).';'r,d> , 

ij  ij 


(3.7) 


(»t,r.d)  .  ’4  <■*»>)  . 

*#3  ifj 


Black  a tat 


(k, black)  (k+1,red)  (k+1,red)  (k+1#red)  (k+1,red) 

‘ij  “  wi+1,j  i-1#j  Wi# j+1  "i,j-1 


(4k)2  , 


(b)  '4  •«•<*)  .  4  „  .  ul.(k.bl«Ok)  _ 


(3.4) 


(c)  .(Wl.bl.ok)  .  w’/i  .buck),  _ 

Bach  stage  can  ba  carried  out  in  parallel#  with  the  red( black)  processors 
working  and  the  black(red)  processors  idle. 

Ibis  idea  of  using  the  red-black  ordering  for  parallel  processors  has  appeared 
several  tinea  in  the  literature  (Heller  [ 1978] ) .  xta  use  on  DAP  was  first  suggested 
by  Hunt  [1974].  (In  Europe#  white-black  chessboards  are  sore  usual  than  red- black 
ones) • 

Zn  Figure  3.5  we  list  a  DAP-FOftTRAN  subroutine  PH0J8OR  for  implementing  the 
heart  of  the  algorithm  (3.7),  (3.8).  The  subroutine  is  provided  with  several  input 
parameters  with  obvious  meanings.  Zn  addition#  two  logical  matrices  are  provided  as 
inputt  the  logical  matrix  MABKMA0K  is  true  at  gridpoints  in  the  interior  of  the 
dam#  and  false  elsewhere i  the  logicsl  matrix  MASK  is  true  at  black  gridpoints  and 
false  at  red  gridpoints.  Finally,  the  values  of  the  real  matrix  w  at  the  boundary 
points  3R  must  be  computed  using  (3.4d)  before  PROJSOR  is  called.  A  full  listing 
of  the  program  is  given  in  Appendix  B. 
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SUBROUTINE  PROJSOR 

COMMON  /RMAT/W(  ,  ) 

COMMON  /RSCA/MAX  DIFF, OMEGA, EPSILON, DAM  WIDTH, DAM  HEIGHT 
COMMON  /I SC A/NUMB  ITERATIONS, NUMB  ROWS, NUMB  COLS 
COMMON  /SUBLMAT/MASK(  ,  ),  MASK  MASK(  ,  ) 

REAL  W,MAX  DIFF , OMEGA, EPSILON , DAM  WIDTH, DAM  HEIGHT 
LOGICAL  MASK,  MASK  MASK 

INTEGER  NUMB  ITERATIONS, NUMB  ROWS, NUMB  COLS 


REAL  Z(  ,  ),  GRID2,  ZMIN  W(  ,  ),SAVEW( 
REAL  ALPHA, BETA 
INTEGER  NUMB  TIMES 
LOGICAL  DONE,  WSIGN(  ,  ) 

EQUIVALENCE  (WSIGN,  W) 

W(WSIGN)  -  0.0 


ALPHA  -  OMEGA  *  .25 
BETA  -  1.0  -  OMEGA 

GRID2  -  (DAM  HEIGHT/NUMB  ROWS)  ‘*  2 
40  SAVE  W  -  W 

NUMB  ITERATIONS  -  NUMB  ITERATIONS  +  1 
DO  45  NUMB  TIMES  -  1,2 

1  MASK (MASK  MASK)  *  .NOT.  MASK 

2  Z  -  W  +  W( -,- ) 

3  Z  -  Z(+,  )  +  Z(  ,+)  -  GRID2 

4  W(MASK)  “  ALPHA  *  Z  +  BETA  *  W 

5  W( WSIGN  .AJO.  MASK  MASK)  -  0.0 
45  CONTINUE 

MAX  DIFF  -  MAX ( ABS ( SAVEW  -  W) ) 


DONE  -  (MAX  DIFF. LE. EPSILON) 
IF  (.NOT.  DONE)  GO  TO  40 

RETURN 


FIG.  3.5.  The  DAP 


,  )  Local  variables. 


Ensure  that  W  is  nonnegative 
everywhere. 

Calculate  the  constants  that  are 
needed  later  on. 


Start  main  loop. 

Save  the  old  value  of  W. 


Reverse  state  of  MASK. 

Calculate  Z  on  only  the  red  (or 
black)  points  as  determined  by 
the  MASK. 

Project 


Find  maximum  difference  between 
old  and  new. 

Check  if  desired  accuracy  is 
attained. 


subroutine  PROJSOR. 
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The  computation  time  for  one  pass  through  the  main  loop  of  the  subroutine 


PROJSOR  is  estimated  in  Figure  3.6,  from  which  it  follows  that  each  PROJSOR 
iteration,  which  requires  two  passes  through  the  loop,  takes  about 
2  x  1135  “  2.27ms.  To  check  this  estimate,  the  average  execution  time  per 
iteration  in  the  subroutine  PROJSOR  was  obtained  by  measuring  (on  a  real  external 
physical  clock)  the  time  required  for  a  large  number  of  iterations  for  the  dam 
problem  with  H  -  24,  h  *  0,  L  ■  16,  and  Ax  -  1.  (This  particular  problem  was 
chosen  because  it  is  a  test  problem  which  has  been  solved  by  many  authors).  The 
measured  time  per  iteration  on  the  pilot  DAP  was  2.2ms,  as  compared  to  the  estimated 
time  of  2.27ms. 


Statement 

Operations 

Time  (ps) 

1 

MASK( MASKMASK )  -  .MOT.  MASK 

1 

logical  mask 

1 

1 

logical  negation 

1 

1 

logical  store 

1 

2 

Z  -  W  +  W(— ,-) 

1 

index  shift  of  two  places 

21 

1 

floating  point  matrix  addition 

160 

3 

Z  -  Z(+,  )  +  Z(  ,+)  -  GRID2 

2 

index  shifts 

30 

1 

floating  point  matrix  addition 

160 

1 

floating  point  matrix 

subtraction 

160 

1 

scalar-matrix  assignment 

15 

4 

W( MASK)  -  ALPHA  *  Z 

1 

floating  point  matrix  addition 

160 

+  BETA  *  W 

2 

floating  point  matrix 

multiplications  by  a  constant 

400 

1 

logical  mask 

1 

5 

W( WSIGN  .AMD.  MASK 

1 

logical  AND 

2 

MASK)  -  0.0 

1 

logical  mask 

1 

1 

scalar-matrix  assignment 

15 

DO 

15  NUMB  TIMES  -  1,2 

7 

1135 

FIG.  3.6.  Estimated  computation  time  for  the  inner  loop  of  PROJSOR. 
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We  conclude  this  section  with  some  comments » 

1.  For  comparison,  the  dam  problem  with  H  -  24,  h  ■  0,  L  -  16,  and  Ax  ■  1  was  also 
solved  on  the  UNI VAC  1180  at  the  University  of  Wisconsin,  using  the  conventional 
ordering  of  gridpoints  and  an  optimizing  compiler  with  single  precision 
arithmetic  (36  bits),  and  the  time  per  iteration  was  found  to  be  5.29ms.  For 
this  problem  the  Pilot  DAP  was  therefore  2.4  times  faster  than  the  UNIVAC  1180. 

It  should  be  noted  that  for  this  problem  only  25  x  17  ■  425  of  the  1024 
DAP  PEs  were  used.  On  a  square  region  the  Pilot  DAP  would  be  six  times  faster 
than  the  UNIVAC  1180. 

2.  In  general,  one  expects  to  be  able  to  predict  DAP  execution  times  to  within 
about  5%,  because  DAP  programs  have  little  overhead  and  spend  almost  all  their 
time  in  computation. 

3.  Since  DAP  floating  point  operations  are  relatively  expensive,  it  is  worthwhile 
optimizing  the  code.  (Readers  who  used  early  computers  which  also  had 
relatively  slow  arithmetic  operations  may  feel  nostalgic) .  An  example  of  such 
optimization  occurs  in  the  subroutine  PROJSOR  (see  Figure  3.5).  The  computation 
(3.7a)  could  have  been  implemented  as: 

Z  -  W(+,  )  +  W(-,  )  +  W(  ,+)  +  W(  ,-)  -  GRID 2 

which  requires  three  additions  and  one  subtraction,  and  takes 

4(15)  +  4(160)  +  15  -  745)as  . 

(shifts)  (additions)  (scalar-matrix  assignment) 

However,  by  sharing  intermediate  xesults  between  PE's,  the  amount  of  arithmetic  can 
be  reduced;  the  implementation  in  PROJSOR  is 

Z  -  W  +  W(  -, - ) 

Z  -  Z<+,  )  +  Z(+,  )  -  GRID2 
which  is  estimated  at  only  546ms. 


It  should  be  noted  that  both  implementations  use  only  half  the  PEs  for 
arithmetic  at  any  one  time.  Larger  grids  or  three-dimensional  problems  (see  Section 
4)  can  use  all  the  PEs  simultaneously. 
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4.  Numerical  solution  of  a  three-dimensional  free  boundary  problem. 

A  three-dimensional  extension  of  the  dam  problem  of  Figure  3*1  was  introduced 
by  Stampacchia  [1974]  (see  also  France  [1974])*  Water  seeps  through  a  porous  dam  in 
a  rectangular  channel  of  width  a  and  height  H •  The  walls  of  the  dam  are  vertical 
but  the  thickness  of  the  dam  is  variable,  so  that  the  dam  occupies  the  region 

X  <0,H)  '  (4.1) 

where  the  horizontal  cross-section  (1^  is  of  the  form 

■  { (x,y)  :  0  <  x  <  a,  ^(x)  <  y  <  V>2(x)}  .  (4.2) 

In  the  specific  problem  considered  here,  is  the  L-shaped  region 

a2  -  ( 0 ,ED )  X  (0,FE)  U  [ED, AF)  x  (0,AB)  ,  (4.3) 

where  the  points  A,  B,  C,  D,  E  and  F  are  as  shown  in  Figure  4.1.  The  upstream 
water  height  is  H  and  the  downstream  water  height  is  h. 

As  shown  by  stampacchia  [1974],  the  problem  can  be  formulated  as  follows: 

Find  u  on  the  region  such  that: 

(a)  -V2u  «  -  [u  +  u  +  u  ]  >  -1,  in  R,  , 

XX  yy  zz  -  3 

(b)  u  £  0.-  in  £l3  ,  (4.4) 

(c)  u(-V2u  +  1)  ■  0,  in  , 
and 


- 

z)2. 

on 

the 

i  <h 

- 

z)2. 

on 

the 

U  SB  g  m  • 

0 

# 

on 

the 

0 

9 

on 

the 

,  a(x,y) , 

on 

the 

upstream  face  AAqFqF  , 
downstream  face  below  water 
downstream  face  above  water 
top  ABCDEF  , 
bottom  A0B0C0D0E0F0 


level  BocoD0E:0E1D1c1B1  , 
level  B^D^EDCB  ,  (4.5) 


and 


u  =  u  «  0,  on  the  sides  ABBnAn  and  EFFnEn  . 
x  n  0  0  0  0 


(4.6) 
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FIG.  4.1.  Flow  through  a  three-dimensional  porous  dam  with 
L-shaped  horizontal  cross- section. 


Here  a(x,y)  Is  the  solution  of  the  two-dimensional  mixed  boundary  value  problem 


(a) 


(b) 


(c) 


a  +  a  =  0,  in  A.B.C-D.EF.  , 
xx  yy  000000' 


2  H2'  on  Vo 


(4.7) 


1  2 

^  j  h  ,  on  b0CqD0E0 


a  =  a 
x  n 


0,  on  A.B.  U  E„F„ 
0  0  0  0 


To  solve  the  problem  (4. 4)- (4. 7)  numerically  we  introduce  a  grid  with 

Ax  *  Ay  «  Az  and  denote  the  approximation  to  u([i  -  2]Ax, [j  -  1]Ay,[k  -  1]Az) 

by  w.j, ,  and  the  approximation  to  a([i  -  2]Ax,[j  -  1]Ay)  by  w. .  i  w.  for 
iJK  13’ 

2  <  i  <  M  -  1  and  1  _<  j  <  N.  As  in  Bruch  [1980]  the  computation  proceeds  in  two 

stages  t 


Stage  1»  The  two-dimensional  problem  (4.7)  is  approximated  by  replacing  the 
differential  equation  (4.7a)  by  the  difference  equations 
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f 


4wij  -  wi+1,j  -  wi-1,j  -  wi,j+1  -  wi,j-1  =  0  • 


(4.8) 


The  Dirichlet  boundary  conditions  (4.7b)  are  satisfied  by  computing  and  storing  the 
values  of  w^j1  =  on  AoF0  B0C0D0E0*  Tlie  Neumann  conditions  (4.7c)  are 

satisfied  by  introducing  two  fictitious  rows  of  gridpoints,  adjacent  to  AqBq  30(1 
EqFq  respectively,  and  requiring  that  the  values  of  w  on  a  fictitious  row  should 
be  equal  to  the  values  of  w  on  the  corresponding  interior  row;  that  is, 

W1  j  =  w3j  and  wM,j  =  wM-2,j'  for  1  <  3  <  N. 

The  resulting  system  of  equations  is  solved  using  a  simple  modification  of  the 

subroutine  PROJSOR  (see  Figure  3.5):  the  term  -GRID2  is  dropped  from  statement 
number  3;  statement  number  5  is  deleted;  and  the  statements 


W(1,  )  =  W( 3,  )  , 

W(M,  )  =  W(M  -  2,  )  , 


(4.9) 


are  inserted  between  statements  number  1  and  2,  so  as  to  make  the  values  at  the 

fictitious  points  equal  to  the  corresponding  interior  values; 

Stage  II ;  The  three-dimensional  problem  (4.4)  is  approximated  by  the  LCP 

(a)  6w  ,  ,  >  w_  .  ,  +  w,  .  ,  +  w.  .  +  w.  .  + 

i,3,k—  i+1,],k  i-1,3,k  1,3+1, k  1,3-1, k 

2 

+  w.  .  .  ,  +  w.  .  ...  -  (Ax)  , 

i,3,k-1  1,3 ,k+1 


(4.10) 


w.  .  .  >  0  , 

- 

wi , j ,k ^Wi , j ,k  Wi+1,j,k  wi-1,j,k  Wi,j+1,k  Wi,j-1,k 

2 

-  w  .  -  w  .  ,  +  (Ax)  ]  =  0  . 

i,3,k+1  i,3,k-1 


The  Dirichlet  boundary  conditions  (4.5)  are  readily  imposed,  while  the  Neumann 
conditions  (4.6)  are  treated  by  introducing  fictitious  sides  parallel  to  the  sides 
ABBqAq  and  EFFqEq,  and  requiring  that  the  values  of  w  on  the  fictitious  sides 
be  equal  to  the  values  of  w  at  the  corresponding  interior  points. 
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To  solve  the  LCP  (4.10)  we  introduce  a  three-dimensional  red-black  partitioning 
of  the  gridpoints,  so  that  each  red( black)  gridpoint  has  six  black(red)  neighbors. 
(It  should  be  noted  that  the  red/black  ordering  on  any  horizontal  plane  is  the  nega¬ 
tion  of  the  red/black  orderings  on  the  adjacent  horizontal  planes.)  As  in  the  two- 
dimensional  problem  treated  in  Section  3,  each  projected  SOR  iteration  can  be  broken 
down  into  two  stages:  a  red  stage  in  which  projected  SOR  is  applied  to  all  the  red 
points  in  the  three-dimensional  w  array,  followed  by  a  similar  black  stage.  In 
detail: 

Red  Stage 


(a) 


(k,red)  (k, black)  (k,black)  (k,black)  .  (k, black) 

ijk  i+1,j,k  i-1,j,k  i,j+1,k  i,j-1,k 


(k, black)  ,  (k, black) 

+  w.  .  +  w.  ... 

1,3, k+1  i,3,k-1 


(Ax)2  , 


(b) 


( k+  V?  ,red) 
w  ^ 
ijk 
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(4.12) 


.  .  (k+1, black)  j-  ( k+  Vo  » black ) •, 

(O  wijk  =  max{0,w.jk  }  . 

To  implement  the  algorithm  (4.11),  (4.12)  it  was  assumed  that  the  dimensions  of 
02  were  such  that  the  gridpoints  on  any  horizontal  cross-section  of  the  dam  could 
be  regarded  as  a  subset  of  a  32  x  32  array.  The  solution  w  was  stored  as  an  array 
of  matrices,  the  matrix  W(  ,,k)  containing  the  values  of  w  on  the  horizontal 
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plane  at  a  height  <k  -  1)Az.  To  control  the  parallel  computation  two  logical 


matrices  were  used:  MASKRB  which  is  true  at  interior  red  gridpoints  in  the  current 
horizontal  cross-section  and  false  otherwise;  and  MASKMASK  which  is  true  at  interior 
points  of  fi2  and  false  otherwise. 

The  algorithm  (4.11),  (4.12)  was  implemented  in  two  ways: 

Implementation  1 : 

During  each  red( black)  stage  the  horizontal  planes  were  updated  in  turn,  and  on 
each  plane  the  red( black)  points  were  updated  in  parallel. 

The  computation  of  z^5^  requires  five  additions  and  one  subtraction.  Given 
an  unlimited  number  of  processors,  n  additions/subtractions  require  log2n  steps, 
so  that  six  additions/ subtractions  require  at  least  three  steps.  By  taking 
advantage  of  idle  PEs,  and  remembering  that,  on  the  DAP,  shift  operations  are  much 
faster  than  arithmetic  operations,  the  DAP-FORTRAN  subroutine  in  Figure  4.2  is  an 
efficient  implementation  of  (4.11),  (4.12)  (compare  Figure  3.5).  A  iull  listing  of 
the  program  is  given  in  Appendix  C. 

The  subroutine  in  Figure  4.2  uses  the  functions  SHS(outh)  and  SHN(orth)  to 
shift  w  instead  of  the  equivalent,  but  slower,  statements  (4.9). 

Implementation  2: 

As  in  the  three-dimensional  magnetohydrodynamic  code  of  Reddaway  [1976]  we  re¬ 
arrange  the  values  of  w.  The  horizontal  planes  are  considered  in  pairs,  and  che 
red  points  on  each  even-numbered  plane  are  exchanged  with  the  corresponding  black 
points  on  the  next  odd-numbered  plane.  As  a  result,  instead  of  having  n  planes, 
each  containing  red  and  black  points  in  a  checkerboard  pattern,  we  have  n/2  planes 
of  red  points  interleaved  with  n/2  planes  of  black  points.  This  makes  it  possible 
to  use  simultaneously  all  interior  PEs  for  arithmetic. 

The  corresponding  subroutine  is  given  in  Figure  4.3,  and  a  full  listing  of  the 
program  is  given  in  Appendix  D.  To  save  time  the  test  for  convergence  is  executed 
only  every  TIMES  iterations,  where  TIMES  is  an  input  parameter. 
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The  subroutine  in  Figure  4.3  assumes  that  there  is  an  even  number  of  planes. 

To  avoid  additional  testing/  it  is  assumed  that  a  copy  of  the  top  plane  is  stored 
above  the  top  plane. 

The  two  implementations  were  run  on  the  problem  with  H  ”  10,  h  ■  0, 

AF  -  FE  -  20,  CD  -  BC  *  10,  which  was  chosen  because  it  had  previously  been  solved 
by  Bruch  [1980].  For  comparison,  the  problem  was  also  solved  on  the  UNXVAC  1180 
using  single  precision  arithmetic  and  optimized  FORTRAN  code.  The  measured  computa¬ 
tion  times  per  projected  SOR  iteration  (including  both  red  and  black  stages)  were: 

Implementation  1:  32ms 

Implementation  2:  16. 0-1 8. 2ms 

(dependent  on  frequency  of  convergence  tests) 

UNI VAC  1180:  34ms 

so  that  implementation  2  on  the  Pilot  DAP  is  about  2  times  faster  than  the  UNIVAC 
1180. 

The  estimated  time  per  SOR  iteration  (implementation  2)  was  found  as  in  Figure 
3.6,  and  was  found  to  lie  between  15.4  and  17.4ms,  depending  upon  the  frequency  of 
convergence  tests. 

For  this  problem  only  383  (i.e.  21  x  23  -  10  x  10)  of  the  1024  PEs  were  used. 
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THE  MAIN  LOOP  -  PROCESS  ALL  THE  Z  PLANES 

SUBROUTINE  MAIN  LOOP 
COMMON  /ISCA/  TOPPLANE ,M 

COMMON  /ISCA/  DAMMAXITERS,  BOTTOMMAXITERS ,  NUMBITERS ,  NUMBOT 
INTEGER  TOPPLANE, M 

INTEGER  DAMMAXITERS,  BOTTOMMAXITERS,  NUMBITERS 
REAL  DAMEPSILON,  BOTTOMEPSILON,  OMEGA,  MAXDIFF 
COMMON  /RSCA/  DAMEPSILON,  BOTTOMEPSILON,  OMEGA,  MAXDIFF 
COMMON  /RMAT/W(  ,,25) 

COMMON  /SUBLMAT/ MASKRB  (  ,  ),  MASKMASK (  ,  ) 

LOGICAL  MASKRB,  MASKMASK 

REAL  SAVEW(  ,  ),  Z  (  ,  ),  Zl(  ,  ) 

REAL  MAXSOFAR,  ALPHA,  BETA,  WIDGRID2,  WIDTHGRID 

INTEGER  NUMBTIMES,  TOPPLANE 

LOGICAL  TEMPMASKf  ,  )  ,  DONE,  WSIGN(  ,  ) 

EQUIVALENCE  (WSIGN,Z) 

ALPHA  =  OMEGA  *  1.0  /  6.0 
BETA  -  1.0  -  OMEGA 

WIDTH  OF  GRID  (I.E.  ONE  UNIT  SQUARE)  IS  SET  TO  1.0 

NUMBITERS  =  0 
WIDTHGRID  =  1.0 

WIDGRID2  =  WIDTHGRID  *  WIDTHGRID 
SAVE  THE  MASKRB  FOR  LATER  RESTORATION 
TEMPMASK  =  MASKRB 

MAXDIFF  IS  THE  MAXIMUM  DIFFERENCE  BETWEEN  SAVEW  AND  W(  , ,K)  AFTER  W(  ,,K) 
HAS  ITS  RED  (OR  BLACK)  VALUES  CHANGE  (FOR  ALL  K) 

MAXDIFF  =0.0 
NUMBITERS  =  NUMBITERS  +  1 
MASKRB  =  TEMPMASK 
DO  30  NUMBTIMES  =  1,2 

ITERATE  FROM  THE  2ND  PLANE  TO  THE  TOP  PLANE 

DO  20  K  =  2,  TOPPLANE 
SAVEW  =  W(  ,,K) 

REVERSE  RED/BLACK  FOR  SUCCESSIVE  PLANES 

MASKRB (MASKMASK)  -  .NOT.  MASKRB 

SUM  THE  SIX  NEIGHBORS 

SAVEW(1,  )  *  SHN ( SAVEW , 2 ) 

SAVEW ( M,  )  =  SHS (SAVEW, 2 ) 

Z  =  SAVEW ( -  ,  - ) 

Z1  -  SAVEW 

Z1 (MASKRB)  •  W{  ,,K  +  1) 

Z ( MASKRB )  *  W(  ,,K  -  1 ) 

Z  »  Z  +  Z1 
Z1  =  Z(  ,+) 


o  o 


Z 1 ( . NOT .MASKRB)  »  -WIDGRID2 
Z  -  Z  +  Z1 
Z  -  Z  +  Z(  +  ,  ) 

C 

C  STORE  THE  AVERAGE  OF  THE  SIX  NEIGHBORS  IN  W  ONLY  IN  THE  RED  (OR  BLACK)  CELLS 
C 

Z  »  ALPHA  *  Z  +  BETA  *  SAVEW 
Z(WSIGN)  =  0.0 
W( MASKRB, K)  =  Z 
C 

C  FIND  THE  MAXIMUM  DIFFERENCE  ON  THIS  PLANE 
C 

MAXSOFAR  =  MAX ( ABS ( SAVEW  -  Z) , MASKRB) 

IF  (MAXSOFAR  .GT.  MAXDIFF)  MAXDIFF  -  MAXSOFAR 
20  CONTINUE 

C 

C  REVERSE  STATE  OF  ORIGINAL  MASKRB  FOR  THE  2ND  PASS  THROUGH  THE  PLANES 
C 

MASKRB (MASKMASK)  =  .NOT.  TEMPMASK 
30  CONTINUE 

DONE  =  (NUMBITERS.GT.DAMMAXITERS). OR. (MAXDIFF. LE.DAMEPSILON) 

IF  (.NOT.  DONE)  GOTO  10 
MASKRB  =  TEMPMASK 
RETURN 
END 


FIG.  4.2.  First  implementation  of  (4.11)  and  (4.12). 
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THE  MAIN  LOOP  -  PROCESS  ALL  THE  Z  PLANES 

SUBROUTINE  MAIN  LOOP 
COMMON  /ISCA/  TOPPLANE , M 

COMMON  /ISCA/  DAMMAXITERS,  BOTTOMMAXITERS ,  NUMB ITERS,  NUMBOT 
INTEGER  TOPPLANE, M 

INTEGER  DAMMAXITERS,  BOTTOMMAXITERS,  NUMBITERS,  TIMES 

REAL  DAMEPSILON,  BOTTOMEPSILON,  OMEGA,  MAXDIFF 

COMMON  /RSCA/  DAMEPSILON,  BOTTOMEPSILON,  OMEGA,  MAXDIFF,  TIMES 

COMMON  /RMAT/W(  ,,25) 

COMMON  /SUBLMAT/MASKRB (  ,  ),  MASKMASK(  ,  ) 

COMMON  /WORK/  Z,WK 
LOGICAL  MASKRB,  MASKMASK 

REAL  SAVEW(  ,  ),  Z{  ,  ),  Z1(,  ),  WKP1(  ,  ) ,  WK(  ,  ),  MAXD(  ,  ) 
REAL  MAXSOFAR,  ALPHA,  BETA,  WIDGRID2,  WIDTHGRID 
INTEGER  NUMBTIME S ,  TOPPLANE 

LOGICAL  TEMPMASK (  ,  ),  DONE,  WSIGN(  ,  ),  TEST,  NOTTEST 
EQUIVALENCE  (WSIGN,Z) ,  (Z,Z1),  (WK,WKP1 ) 

ALPHA  =  OMEGA  *  1.0  /  6.0 
BETA  =1.0-  OMEGA 
C 

C  WIDTH  OF  GRID  (I.E.  ONE  UNIT  SQUARE)  IS  SET  TO  1.0 
C 

NUMBITERS  =  0 
WIDTHGRID  =  1.0 

WIDGRID2  =  WIDTHGRID  *  WIDTHGRID 
C 
C 

1  TEST  =  TIMES  .EQ.  1 
NOTTEST  =  .NOT.  TEST 
ITIMES  *  TIMES 

MAXD  =  0.0 
C 
C 

C  ALTER  ALL  THE  ODD  NUMBERED  PLANES: 

2  KM2  =  1 

DO  20  K  =  2, TOPPLANE, 2 
SAVEW  =  W(  ,,K  +  1) 

WK  =  W(  ,,K) 

WK( 1 ,  )  =  SHN ( WK, 2 ) 

WK(M,  )  =  SHS( WK,2) 

Z  =  WK  +  WK(  - ,  -  ) 

Z  =  (Z(+,  )  +  Z(  ,+)  +  WK  +  MERGE (W(  ,,KM2),W(  ,,K  +  2), MASKRB) 
-  WIDGRID2)  *  ALPHA  +  SAVEW  *  BETA 
Z(WSIGN)  =0.0 
W( MASKMASK, K  +  1 )  =  Z 
IF  (NOTTEST)  GOTO  20 
Z1  =  ABS( SAVEW  -  Z) 

MAXD(ZI.GT.MAXD)  =  Z1 
20  KM2  =  K 

C 
C 
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C  ALTER  ALL  TBI  EVEN  NUMBERED  PLANX81 
DO  21  K  ■  2 , TOPPLANE , 2 
SAVEN  -  W(  ,,K) 

WKP1  -  W(  ,,K  +  1) 

HKP1 ( 1 ,  )  -  SHN(WKP1,2) 

WKP1 (M,  )  -  SBS ( MRP  1,2) 

Z  -  WKP1  ♦  WKPl 

Z  -  (Z(  +  ,  )  +  Z(  ,+  )  +  WKPl  +  MERGE  (W(  ,,K  +  3),W(  ,,X  - 
-  WZDGRID2 )  *  ALPHA  +  SAVEW  *  BETA 
Z(WSXGN)  -  0.0 
W(MASKMASK,K)  -  Z 
IP  (NOTTS ST)  GOTO  21 
Z1  -  ABS(SAVEW-Z) 

MAXD(Z1 .GT.MAXD)  -  Zl 
21  CONTINUE 

C 
C 

ITIMES  -  ITIMES  -  1 
IP  ( ITIMES. GT. 1 )  GOTO  2 
IP  ( ITIMES. EQ.O )  GOTO  3 
TEST  -  .TRUE. 

NOTTEST  -  .FALSE. 

GOTO  2 
C 

3  NUMB  ITERS  -  NUMB  ITERS  +  TIMES 
MAXDIPP  -  MAX ( MAXD , MASKMASK ) 

IP  (MAXDIPP. LT.DAMEPSILON)  GOTO  4 
C 

IP  (NUMBITERS.LT. DAMMAXITERS)  GOTO  1 

4  RETURN 
END 


FIG.  4.3.  Second  implementation  of  (4.11)  and  (4.12). 
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(a)  For  purposes  of  comparison  we  have  used  previously  published  problems  but  they 
have  dimensions  which  do  not  match  the  DAP  array  closely.  In  many  practical 
problems  the  resolution  would  be  tailored  to  the  DAP  dimensions  to  achieve 
higher  performance. 

(b)  The  programs  presented  are  readily  extensible  to  larger  problems  on 
correspondingly  larger  DAPs  such  as  the  production  64  x  64;  it  is  only 
necessary  to  change  the  boundaries.  The  time  to  process  one  plane  would  be 
unchanged. 

(c)  Performance  on  small  three-dimensional  problems  can  be  improved  by  mapping 
several  problem  planes  onto  one  DAP  matrix. 

(d)  Problems  with  large  "horizontal"  dimensions  can  be  mapped  with  each  PE  holding 
a  small  neighborhood  group  of  points.  Performance  improves  because  each  PE 
holds  both  black  and  red  points  (Hunt  {1979]). 

(e)  Very  large  problems  cannot  be  held  entirely  within  DAP  store.  For  example  with 
four  times  as  many  points  in  a  horizontal  plane  as  there  are  PEs  the  limit  is 
about  26  planes  with  4K  bits  per  PE  or  about  122  planes  with  16K  bits  per  PE. 
With  backing  store  the  transfer  rates  with  N  active  problem  planes  in  the  DAP 
can  be  minimized  by  advancing  each  plane  (N  -  2)/2  iterations  per  backing 
store  fetch.  Hence  it  should  be  possible  to  achieve  a  balance  between  input- 
output  and  processing  times  (Reddaway  [1976]). 

(f)  Problems  of  this  type  offer  possibilities  for  using  fixed  point  arithmetic 
(with  suitable  scaling)  and  using  low  precision  for  computing  the  iterative 
corrections.  This  is  much  faster  than  floating  point  work  and  performance 


improvements  as  large  as  a  factor  of  10  are  predicted  without  loss  of  accuracy 
in  the  final  solution. 
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6.  Conclusions 

He  have  demonstrated  that  two-  and  three- dimensional  linear  complementarity 
problems  can  be  solved  on  DAP  with  high  performance  and  easy  progr earning  using  a 
version  of  projected  SOR.  There  is  scope  for  even  higher  performance  and  for 
tackling  a  wide  range  of  problem  sizes* 
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APPENDIX  A:  The  Pilot  Host- DAP  Interface. 

In  the  Pilot  DAP  system  the  store  of  the  DAP  is  not  an  integral  part  of  the 
host's  store  as  with  the  production  DAP's.  It  is  therefore  necessary  to  explicitly 
move  data  between  the  host  and  DAP  and  this  is  achieved  by  using  standard  host 
FORTRAN  subroutines.  The  subroutine  names  begin  with  DAPTO  or  DAPFROM  depending  on 
whether  they  move  data  into  or  out  of  the  DAP.  The  remaining  letters  of  the  name 
indicate  the  type  (integer  or  real  denoted  by  I  or  E)  and  rank  (scalar,  vector  or 
matrix  denoted  by  S,  V,  or  M)  of  the  variable  transferred.  Parameters  of  DAPTO  and 
DAPFROM  give  the  name  of  the  host  program  variable  and  the  location  within  the  DAP 
in  terms  of  the  name  of  the  common  area  and  the  offset  from  the  start  of  this  area. 

Initiation  of  DAP  processing  is  also  less  direct  on  the  Pilot  system  with  DAP- 
FORTRAN  subroutines  being  called  via  the  standard  host  FORTRAN  subroutine  DAPGO.  A 
statement  of  the  form: 

CALL  DAPGO ( ' DAPSUB '  ,  N ) 

will  suspend  execution  of  the  host  FORTRAN  and  transfer  control  to  the  DAP-FORTRAN 
subroutine  DAPSUB.  Execution  of  the  host  FORTRAN  is  resumed  after  DAPSUB  and  any 
further  levels  of  DAP-FORTRAN  subroutines  have  been  executed.  The  parameter  N 
gives  the  maximum  number  of  seconds  allowed  for  DAP  processing. 


-35- 


n  o 


APPENDIX  B:  The  two-dimensional  dam  problem. 


WASTER  ECALpAL 
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C  RESEARCH  AND  ADVANCED  DEVELOPMENT  CENTRE  IN  STEyENAGE 

c  England  ijune  i/«zo,  1979), 
c 

C  E APL A  I  NAT  I  ON  OF  VARIABLES 

C  w  -  T  IF  VALUES  AT  THE  GRIOPOINTS 

t  Z  -  THE  '*FW  VALUES  AT  THE  GRIOPOINTS 

C  Mas*  .  USED  TO  IMPLEMENT  THE  PEO/BLACk  SCHEME 

C  MSS*.  •IAS"'  -  UsEU  IN  SETTING  Up  THE  MASK*  and  IN  SWITCHING  STATES 

C  ACC  EppjS  -  DIFFERENCES  BETWEEN  THE  OLD  AND  THE  NEW  VALUES 

C 

C  0 *£ SA  -  | H£  OveHRELAxAtIUN  PARAMETER 

c  epsilon  «  the  desired  accuracy 

C  DONE  -  h  TE..1PORART  LOGICAL  VARIABLE 
C  NJMB  IIc^mTI DNS  -  the  NUMBER  UF  ITEHATIUNS 
C  MAX  ITERATIONS  «  THE  MAXIMUM  NUMBER  OF  ITERAUUNS 
C  N-.MA  RO"5>  -  The  NUMBER  OF  HOWS  IN  THE  GRID 

C  N.IMJ  COLS  -  The  NUMBER  OF  COLUMNS  IN  THE  sr ID 

C  Mpln  G-'l"  -  WIDTH  OF  ONE  UNIT  IN  ThE  GKId 
C  HEIGHT  GRID  i  HEIGHT  OF  UNE  UNIT  IN  THE  GRID 
C  JAM  WIDTH  -  SEL‘C  ExPLANAT ORY 
C  DAM  HEIGHT  -  SELF  EXPLANATORY 

C  x  GRIOPOINTS  -  THE  NUMbfcH  OF  GR1UP0INT  IN  THE  X  UiREClIjN 

C  r  f  J  R I  p  0 1 N  T  s  -  THE  NUMBER  of  griopoint  in  the  r  DIRECTION 

C  'iTTE  J  IT  IS  ASSUMED  THAT  HEIGHT  GRID  ■  WIDTH  GRID 
C 

C  uffCLA  RA  fl  ,)N  JF  VARIABLES 
C-M.'.Ow  /p.1AT/w<  32*32  ) 

CD.r'J.  /rSCA/MAX  DIFF,  OMEGA,  EPSILON,  DAM  W  iU  T  H ,  DAM  HEIGHT 
COMMON  /ISCA/MAX  ITERATIONS,  NUMB  ITERATIONS,  NUMB  ROWS,  NUMp  CJLS 
1  <  ^R I DpQ 1  NTS , Y  GRIOPOINTS 

wEAL  •<  t  OMEGA ,  EPSILON,  0  AM  WIDTH,  DAM  HEIGH!  ,  MAX  U1FF 
1MT£„c«  max  ITERATIONS,  NUMB  ITERATIONS,  numb  RUwS,  nJMb  COLS, 
i  *  M.viuP'jiNTS,  Y  GPIDPOIMTS 
PAJSF  *9 

UP.,)  Ii  UMc  NSIJN$  OF  THE  OAM,  AM)  THE  MaXIMjm  NUMBER  UF  DERATIONS 
1  -EA  K?,  10j,ENO»2D)  DAM  WIDTH,  UAM  HEIGHI,  MAX  IlcRAUDNS 

i  00  r  .jc<  'm  T  \  2F  1  w  ,  5 ,  15) 

iIJTh  ,GT,  1,0)  .AND,  (0AM  HEIGHT  ,GI,  l,J)>  GOTO  5 
WO  >^ITF*°*ii  •)  DAM  WIDTH,  UAM  nfciGHT 
t  1 0  f::h  iaM  pr'l  *X  *  *F10,5,8H  NY  ■  fFlQ.S) 

r>  u  ;  t  v 


-36- 


o  n  »-  •*  ui  n  n 


READ  IM  I  HE  Number  OF  RQN$  AND  COLUMNS  in  the  OHIO#  *NO  EFsILON 
»EA0(5»120)  NUMB  CQlSi  NUMB  RUNS,  EPSILON 
20  F0RMAT12I5,  F10,5) 

wRITE<®,lJO)  DAM  HEIGHT,  OAM  WIDTH,  EPSILON*  NUMB  RUwS, 
t  NUMB  COLS,  MAX  ITERATIONS 

30  FORMAT*  17HII  \PU  T  PARAMETERS,  /  14H  DAM  HEIGHT  ■  #E13,P» 
j  l UH  DAM  WIDTH  .  *  E  j  3 . 6 1 1 1 H  EPSlLUN  ■  ,E13,6,/ 

2  I  OH  NUMBER  OF  ROWS  .  ,I5»21H  NUMBER  UF  COLUMNS  ■  *19* 

i  MAXIMUM  ITERATIONS  •, 19) 

concert  *nd  transfer  data  jo  oap 
call  DAPTO  ES(EPS!L3N,4HRSCA,2) 

CALu  DAPTD  ESIOAM  H I u  TH  *  4MRSC A  *  3  ) 

CALL  OAPTO  ESIOAM  HEIGHT, 4HRSCA, 4) 

CALL  UMPTO  Is<MAX  ITERATIONS, 4HISCA, 0  ) 

CALl  OAPTO  IS«NjMB  ROWS, 4HISCA1 2) 

CALL  OAPTO  IS<NUM3  COLS , 4H I  SC A , 3  ) 

c 

c  call  main  dap  furtran  SuBROutIn£ 

CALL  UMPGd( 7HLAPLACE, 10  ) 

C 

c  concert  and  transfer  DAP  data  BACK  noo 
CALL  OApfRDM  EM(w,4HRMAT,0) 

CALl  uAPFHUM  tS( OMEGA, 4HRSCA, 1 > 

CALL  OAPFRDM  ES(MAX  OIFF, 4HRSCA, o> 

CALL  OAPFROM  1 S ( NUMB  I  TER AT  I ONS , 4HI SCA , 1 ) 

CALL  OAPFROM  1 S ( Y  GRIDKOINTS,  4HISCA,  5) 

CALL  OAPFROM  1 Sf  X  GR10P0INTS,  4HISCA,  4) 

C 

c  write  oui  the  scalers 

^RlTE(O,l40)  NUMB  ITERATIONS,  OMEGA,  MA*  OIFF 
140  FORMAT! •  «,///,»  NUMBER  OF  ITERATIONS  *  *,I9,»  OMfeG*  a  ',£13,9, 
l  ,  MAXIMUM  DIFFERENCE  ■  *,E13,B) 

C 


C 

If. rite  Oul  w 

J  ■  Y  mr  1 OPO I NTS 

1 J 

vRITE)°,J50)  l  w  (  1 , U  ) , 

1*1, X  GHlOPOlNTS) 

15: 

FORMAT!'  'iREU.b) 

J  «  J  "  1 

IF  C  )  , G€ ,  1)  Gq I 0  10 

v>ITEI°,17^) 

HO 

FORMAT! •  «,///) 

,)0  3-J  *  =  I,  I  GRIDPUINTS 

^RITE(O,10u)  !*<t,J), 

0  ■  1,  X  QRIOPOINTS) 

i  1  ■- 

r-H-IAM*  '  ,  QF  1  3 , 6  ) 

30 

C  •'  <T  J  -UF 
'.'TO  1 

f*  • 


-37- 


j 


n  n  n  r>  n  n  o  ooooo  no  non  on  no 


C 

C 


C 

C 


SJBROUUNE  LAPLACE 

SUBROUTINE  TO  C*R«Y  OUT  THE  CALCULATIONS  UN  T**  OAP 
COMUON  /RMAT/W{ , ) 

COMMON  /RSCA/mAK  DIFF,  OMEGA,  EPSILON,  UAM  WIDTH,  0 AM  HEIGHT 
COMMON  /ISCA/MAX  ITERATIONS,  NUMB  1TERAUONS,  NUIM*  ROWS*  NUM* 

common  /isca/x  griopoints,  y  griopoints 

rEA|.  mi  MAX  DlFF,  OMEGA,  EPSILON,  0AM  WIDTH,  OAM  HEIGHT 
INTEGER  MAX  ITERATIONS,  NUMB  ITERATIONS,  NUMB  RUNS,  NUMB  COLS 
INTEGER  X  tiR I OP  0 1  NT  S ,  Y  GRlOPOiNTS 

SET  OP  OMp  FORTRAN  COMMON  AREAS 

COMMON  / SO 9RSCA/ HEIGHT  GRID,  WIDTH  GRID 
COMMON  /SUtfLMAT/MAS<<  ,  ),  MASK  MASM,) 
real  HfcIGHT  GRID,  WIDTH  GRID 
logical.  MASK,  MASK  MASK 

INITIlU6  all  VARIABLES 

call  i*it  others 
call  i"JT  mask 

CALL  i'^IT  w 


COLS 


PERFORM  I  he 
CALL  "MIN 


ACTUAL  CALCULATIONS 
LOOP 


WE  HAVE  hqw  EITHER  ACHEIVEO  THE  DESIRED  ACCURACY  U,E,,  MAX  ERR 
<■  EPSILUN)  OR  WE  HAVE  ITERATED  MORE  THAN  MAX  ITERATIONS  TIMES, 
RETURN 

SUBROUTINE  INI T  MASK 

SUBROUTINE  to  INITIALIZE  the  mask 

COMMON  /ISCA/MAA  ITEKAIIONS,  numb  ITEHAI ions, NUMB  rows, numb  CULP 
COMMON  /ISCA/X  gridpoints,  Y  GHIOPOINTS 
common  /SUblmat/mas<( , > ,  mask  Mask ( , i 

INTEGEh  MAX  ITErAIIONS,  numb  ITERATIONS,  NUMB  ROWS,  NUMB  COLS 
INTEGEH  X  GRIUPOINTS,  T  GRIOPUiNTS 
logical  mask,  mask  mask 


MASK  MASK  IS  A  LOGICAL  MASk  USED  IN  SETTING  UP  ANU  IN  CHANGING  THE 
STATE  of  T Hh  Mask,  MASK  MASMI,J)  IS  FALSt  wHENE  I  »  1  OR  I  > 

NUMB  ROWP;  AND  WHERE  J  ■  1  UR  wHfcRE  J  >  NOMH  COLS,  IT  IS  TRUE 
ELSEHHERt, 

MASK  mask  ■  .NOT,  UlTCINUMA  HOWS*  ,0R,  ALTR1N0MB  CUlSJ) 
mask  MMSK< 1 , )  ■  .FALSE, 
mask  ih$K(, 1)  a  ,F*LSE, 

MASK  IS  M’»CH  THE  SAME  AS  MASK  MASK,  ExCEPt  THAI  WHERE  MASK  MASK 

is  true,  mas*  alternates  between  rwut  *nd  false,  trending 
On  the  FOLLOWING! 

if  MASK  MASK ( I  ,  U )  >  TRUE  ANU  l  *  J  IS  fcVEN  I  HE  N  MASK(J,J) 

■  TRUE,  ELSE,  WASK(I,U)  ■  FaCSE, 

MASK  a  ALTC(l)  ,LEQ,  ALTRU) 

MASK( .WOT.MASX  MASK)  a  .FALSE, 

TRACE  127  I  MERGE! 1, u, MASK  MASK!) 

TRACE  *27  IMEHGEI I,U,M«SK)  ) 

return 
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SJBROU I  INK  InIT  OTHERS 


c  Subroutine  to  initialize  the  other  variables 

common  /RSCA/MAX  diff, omega, epsilon, oam  WIDTH, OAM  height 

COMMON  /ISCA/MAX  ITERATIONS, numb  ITERATIONS, numb  hows, numb  COLS 

COMMON  /ISCA/X  GRIDPOINTS,  V  GRIOPOINTS 

COMMON  /SUBRSCA/HEIShT  GR1U,  WIDTH  GRID 

REAL  MAX  OiFF,  OMEGA,  EPSILON,  DAM  WlOTH,  DAM  HEIGHT 

REAL  HftIGHT  OHIO,  WIDTH  GRID 

INTEGEH  X  GRIOPOINTS,  Y  GRIDPOINTS,  MAX  ITERATIONS,  NUMB  ITERATIONS 
INTEGER  NUM§  ROWS,  NUMB  COLS 
C 

c  for  the  »est  run, 

C  dam  HEIGHT  •  2*,0 

C  OAM  WIDTH  •  16,0 

X  GRIOPOINTS  ■  NUMB  COLS  ♦  1 

Y  GRIOPOJNIS  ■  NUMB  ROWS  ♦  1 

HEIGHT  GRID  a  DAM  HEIGHT  /  NUMB  ROWS 

WIDTH  v»RId  ■  uAM  wIdTH  /  NUMb  COLS 

OMEGA  ■  1,8 

NJM9  I  IERATIONS  *  0 

RETURN  _ _ 

SUBROUI INE  INIT  W 

COMMON  /RSCA/MAX  0 I FF , UMEGA , EPS I LON , DAM  WIDTH, DAM  HEIGHT 
REAL  max  DIFF, OMEGA, EPSILON, OAM  WIDTH, DAM  HEIGHT 

COMMON  /ISCA/MAX  ITERATIONS, NUMq  I Te«AT iUNS, NDMB  ROWS , NOM§  gULS 
COMMON  /IsCA/X  GrIDPOINts,  Y  GRIDPOINTS 
COMMON  /SUBRSCA/HEIGHT  GRID,  WIDTH  GRIO 
REAL  HtIGHT  GHID,  WIDTH  GRID 

INTEGE*  MAX  lTERATlONSi  NUMb  I  TgRA  T I  DNS,  NUMa  ROMS,  NUMjj  cOLS 

INTEGER  x  GrIUPOINTSi  t  GrIOPOINTS 

common  /RMAT/W( ,  ) 

real  «•  T£MPV(  ),  TEMPS 

REAL  INOEXi) 

C 

C  INITIALISE  THE  W  MATRIX,  I.E,,  SET  THE  CUNStANj  BOUNDARY  VALUES 
C  And  THE  INITIAL  NQN»BOuNUAR Y  VALUES, 

C 

c  SET  UP  a  VECTOR  CONTAINING  t«e  INDICES 

DO  30  *  ■  1,3* 

INDEX(I)  ■  I  •  1,0 
30  CONTINUE 
c 

C  SET  All  UF  M  ■  ZERO 
W  ■  O.u 

c 

C  CALCULAf6  THD  TfeMPO«AHY  CONSTANlS 

TEMPS  ■  DAM  HEIGHT  •  QAm  HEIGHT  *  0,5 
TEMPv  *  INDEX  /  NUMB  HOWS 

Will)  *  TEMPS  *  (1,0  ■»  INDEX  *  WIDTH  GRlu  /  DAM  wiDTH) 

HU,)  ■  TEMPS  *  (1,0  a  TEMPV)  *  (1,0  *  IEmPV) 

TRACc  127  <W) 

TRACE  127  (TEMPS,  DAM  height,  *10TH  GH 1 o ) 

TRACE  127  (0AM  WIDTH,  height  GRIO) 

TRACE  127  (TFMPV) 
trace  *27  Unde*) 
return 


on  non  no  4>  o  o  o  n  o  oo  on 


SJBPOJI INE  MAIN  LOOP 


c  Subroutine  to  carry  out  the  calculations 

COMMON  /RMAT /«(,) 

COMMON  /RSCA/MAX  OIFF, OMEGA, EPSILON, 0AM  h10TH,0AM  HEIGHT 
COMMON  /IScA/MAK  XTgRAf I0NS,NUMb  l TeRAT IONS, NUMB  KO*$,NUM§  cOL* 
COMMON  /ISCA/*  GrIDPOINTS#  v  gridpojnts 

COMMON  /SO«LMAT/MASKI ,  ),  MASK  MASK ( ,  ) 

COMMON  /SUBRSCA/ HEIGHT  GRIO,  WIDTH  GRID 

REAL  HfMAX  DIPFiOMEGA, EPSILON, 0AM  HIOTHfOAM  HEIGHT 

REAL  MtIGHT  GRID,  WIOTH  GRID 

logical  mask,  mask  MASK 

integen  x  griopoints,  r  griopuints 

INTEGER  MAX  ITERATIONS,  NUMB  ITERATIONS,  numb  ROMS,  NUMB  COLS 
LOCAL  VAHiABi.ES 

REAL  zl,),  WIU  GRIO  2,  z  MIN  W(,l,  SAVE  Mi,) 

INTEGER  NUMB  TIMES 
LOGICAL  DUNE,  M  SIGN( , j 
EQUIVALENCE  (H  SIGN,  HI 

ENSURE  THAT  N  is  >■  ZERO  EVERYWHERE 

W(WSISN,  >  0.0 

CALCULATE  the  CONSTANTS  IHAT  ARE  NEEDED  LATER  UN 
ALPHA  ■  OMEGA  *  ,25 
BETA  ■  1,0  .  OMEGA 

W ID  GRID  2  ■  WIDTH  GRID  *  WIDTH  GRID 

START  MAIN  LOOP 
SAVE  THE  OLD  VALUE  OF  W 
0  save  W  *  W 

NJM*  DERATIONS  a  NJMB  ITERATIONS  ♦  1 
DO  45  NUMB  TIMES  a  1,2 

reverse  *TATE  qP  mask 

mASK( M«SK  MAS*)  •  , NOT ,  MASK 

CALCULAT6  z  UN  UNLY  THE  RED  <0R  BLACK)  PUlNyS  AS  DETERMINED  BY 

the  mask 

Z  ■  N  *  *<•»-•) 

7  ■  7<*»  >  ♦  7<  ,♦>  »  WIU  GRIP  2 
w(mask)  a  alpha  *  l  ♦  beta  *  w 

special  part»  ensure  that  m  is  >■  o 

W(W  Sl^N  , ANp ,  MASK  MASK)  a  0,0 
45  CONTINUE 
C 

C  FInO  MAXIMUM  DIFFERENCE  BETWEEN  Old  ANO  Nkw  a 
VAX  ulfF  «  M  A  X ( A  B  S ( S  A  V  t  H  m  Ml) 

C 

C  CHECK  jO  SEE  IF  DESIRED  ACCURACY  IS  ATTAINED,  UH  IF  NUMBER  OF 
C  ITERATIONS  has  EXCEEDED  LIMIT 

OONEalMAX  DIFF,LE,EPSILON),OR,(NUMB  ITERAT IUNS,G1 ,MAX  I TERA  T I  ONE ) 
IF  (.NUT,  UDNfc)  GO  TO  4o 

C 

TRACE  125  <*| 
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TRACE  *25  (VIA*  DIFF) 

KIUOCV  ,i. 

!F<MUM0  I  T£RaTIONS|ST|MAX  I TEKAT  IONS  JNU*B  i  TEKATiaNf«MA*  ITERATION 

all  oone  main  loop 

RETURN 


ooo 


APPENDIX  C:  The  three-dimensional  dam  problem  -  implementation  1. 


MASTER  ECAWPAL 


10 


COMMON  /rmAT/W<32#32) 

COMMON  /{SC*/  ZPLANES#  XPOINTS#  YPOlNTS.  XFPOINTS#  XBPOINTS 
COMMON  /ISCA/  YRPOINTS,  YLPOINTS,  DAMHEIGHT,  OAMFACEi  OAMRSIOE 
COMMON  /ISCA/  DAMLBACK#  OAMRBACK,  OAMLFSIOE,  OAMLBSIOE 
COMMON  /ISCA/  DAMMAXITERS#  BOTTOMMAX ITERS,  NUMBITERS#  NUMBOT 
1NTE66R  ZPLANES#  XPOINTS.  YPOlNTS#  XFPOINTB#  X|POINT$ 

INTEGER  YRPOINTS#  YLPOINTS^  DAMME JGHT #  DAMPACE#  OAMRSIOE 
INTEGER  DAMLBACK,  DAMRBACK,  DAMLFSIOEi  OAMLBSIOE 
INTEGER  DAMMAXITERS*  BOTTOMMaX I TERS ,  NUMBlTeRS 
PEAL  OAMEPSILON,  80TTQMEPS1L0N#  OMEGA#  MAXOIFF 
COMMON  /RSC A/  OAMEPSILON#  BOTTOMEPS ILON#  OMEGA#  MAXOIFF 
INTEGER  VU32) 

REAL  V R ( 321 

EQUIVALENCE  ( V I ( I),7PLANES) , ( VR<  i ) (OAMEPSILON) 

PAUSE  99 

REA0(2,1C0)  DAMHEIGHT.  OAMLFSIOE.  OaMLBSXDE,  DAMLBACK,  DAMRBACK 
IF(DAMHEIGHT  .LE.  Oj  pause  00 


REAO(2,ltO)  MAXlTERSt  EPSILON 
OAMEPSILON  -  EPSILON 
BQTTOMEPSILON  ■  EPSILON 
DAMMAXITERS  ■  MAXlTERS 
BOTTOMMAXITERS  •  MAXlTERS 


NUMBITERS  ■  1 

OMEGA  i  1.8  _ 

ZPLANES  ■  DAMHEIGHT  ♦  1 
DAMFACE  ■  DAMLBACK  ♦  DAMRBACK 
OAMRSIOE  *  DAMLFSIDE  ♦  OAMLBSIOE 
XPCINTS  «  OAMRSIOE  ♦  1 
XFPOINTS  ■  OAMLFSIOE  ♦  1 
XBPOINTS  »  OAMLBSIOE 
YRPOINTS  «  DAMRBACK  ♦  2 
YLPOINTS  *  DAMLBACK+1 
YPOlNTS  *  YRPOINTS  ♦  YLPOINTS 

WRITEC6.120)  DAMHEIGHT,  OAMLFSIOE.  DAMLBSIDE,  DAMLBACK,  DAMRBACK 
CALL  DAPTO  IV(VI,4HISCA,0) 

CALL  DAPTO  E$< OAMEPSILON, 4HRSCA , 0 > 

CALL  DAPTO  ES( BOTTOMEPSILON,4HRSCA, 1 ) 

CALL  OAPTO  ES( OMEGA ,4HRSCA , 2  ) 


The  time  d‘*t  perioo  is  set  To  boo  seconds 


call  DAP60(7HLAPLACE,  600) 

CALL  oapfrom  FV< Vfi,4MR$CA,0) 
call  DApFROM  I  V < V 1 , 4H I  SC A , 0 ) 

HR  I TE( 6, 1 30 )  NUMBITERS,  NUMBOT# 


MAXOIFF 
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r,.,  ^ _ , .  _  .  ......  . 

C  WRIT?  OUT  THE  w  MATRIX,  FROM  TOP  TO  BOTTOM  PLANE 

nx«  z planes  &  "  . 

M.  NJ  9-  VP 01  NTS  ...  ^ 

nj  ■  xpoinTS 

00. 40  K  »  I, NR  _ 

_  WRITE  (6,150)  R 

CALL  OAPPROM  EM(  N,4HRMAT,32*(K*jl) 

_  00  40  I  •  1 f N1 

iiaiTg.  i 

J0  _  CONTINUE 

goto  io  ..  .  __ 

00  50  J  ■  1 i N J 

WRITE  (0,160)  J  ...i _ i~^=y  i _ . 

DO  50  KK  >  1,NK 

K  m  NKmKKti  . 

CALL  OAPPROM  EM(W,4HRMAT ,32*(R»i ) ) 

'L-L.  -^JiRlTft  (0,4*01  <Wd .H»I*t i Nl.t_^ 

50  CONTINUE 

DO  6Q  I  ■  1,NI 
WRITE  (6,170)  I 
00  60  KKal  ,  NK 
K«NK-KR*t 

1-  .  _£ALL  OAPPROM  EM( W,4MRMAT f32*(K*l)) 

WRI’P  (6,140)  ( W( 1 , J) , J*1 , NJ ) 
to,.  CONTINUE 

GOTO  10 

100  FORMAT! 510) 

110  FORMAT! 10,  EO.O) 

120  FORMAT!  9M1HEIGHY  *,I5, 

1 22H  DAM  LEFT  FRONT  SIDE  ■,I5,20HDAM  LEFT  RACK  SIDE  a, 15 
_  212M  LEFT  BACK  ■, 15, 12MRI0HT  BACK  *,I5#//) 

130  FORMAT ( 23H  NUMBER  OF  ITERATIONS  s,l5, 

19H  NUMBOT  a,I5,20HMAXlMUM  DIFFERENCE  •iEIS.O) 

140  FORMAT! 1M  ,9E13.6) 

150  FORMAT! ///I 1H  FOR  PLANE  ,13) 

170  FORMAT! ///11H  FOR  SIDE  ,13) 

160  FORMAT! /// 1 1H  FOR  FACE  ,13) 

END 
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J+.  1  -  ' 


i 


C  START  Of  0 AP^ORTR A N  SECTION  _ _  1.  .  _ ' 

C  SUBROUTINE  TO  CALL  OTHER  SMSROUTINES  .  I  '  ' 

SUBROUTINE  LAPLACE  7"7~  1T~ 

common  /isca/  {planes,  xpoints*  ypoints,  xfpointe,  xbpoints 

COMMON  /ISCA/  YRPOINTS,  YLPOINTS,  DaMHEIGHT,  DAMFACE,  OAMRSIOE 
COMMON  /1$Ca/  DAMLBACK,  DAMRBACK,  damlfside,  oamlbsioe 
COMMON  /ISCA/  OAMMAXlTERSr  BOTTOMMAX ITERS,  NUMB1TERS,  NUMBOT 
'  INTEGER  ZPLaNES,  XPOINTS,  YPOINTS,  XFPOINTS,  XBPOINTS 
INTEGER  YRPOINTS,  YLPOINTS,  JJAMHEIGhT,  DAMFACE,  OAMRSIDE 
INTEGER  OAMLBACK,  OAMRBACA,  OAMLFSlDE,  DAMLBSIOS 
INTEGER  DAMMAXITERS,  BOTTOMMAX  ITERS,  NUMBITERS 
REAL  OAMEPSILONi  BOTTOMEPSILON,  OMEGA,  MAXOlF? 

COMMON  /rsca/  damepsilon,  bottomepsilon,  omega,  maxoiff 
COMMON  /RMAT/WC , ,25) 

COMMON  /SUBLMAT/MASKRB( ,  ) ,  MASKMASK! ,  ) 

LOGICAL  MASKRB,  MASKMASK 
CALL  INIT  MASK 
CALL  INIT  W 

CALL  INIT  BOTTOM  PLANE 
CALL  MAIN  LOOP 
RETURN 
ENO 

initialize  the  masks 

SUBROUTINE  injt  mask 

COMMON  /ISCA/  ZPLANES,  XPOINTS,  YPOfNTS,  XFPOINTS,  XBPOINTS 
COMMON  /ISCA/  YRPOINTS,  YLPOINTS,  OAMHEIGHT,  OAMFACE,  OAMRSIOE 
COMMON  /ISCA/  DAMLBACK,  OAMRBACK,  DaMLFSIOE,  OAMLBSIOE 
COMMON  /ISCA/  DAMMAXITERS,  BOTTOMMAX I TERS ,  NUMBITERS,  NUMBOT 
INTEGER  ZPLANES,  XPOINTS,  YPOINTS,  XFPOINTS,  XBPOINTS 
INTEGER  YRPOINTS,  YLPOINTS,  OAMHEIGHT*  DAMFACE,  OAMRSIOE 
INTEGER  DAMLBACK,  OAMRBACK,  DAMLFSIOE,  OAMLBSIOE 
INTEGER  DAMMAXITERS,  BQTTOMMAXlTERS,  NUMBITERS 
peal  DAMEPSILON,  BOTTOMEPSILON,  OM^SA,  MAXOIFF 
COMMON  /RSCA/  DAMEPSILON,  BOTTOMEPSILON,  OMEGA,  MAXOIFF 
COMMON  /RMAT/HI , ,25) 

COMMON  /sublmat/maskrbi ,  MASKMASK(,) 

LOGICAL  MASKRB,  MASKMASK 
LOGICAL  T (  ,  ) 

MASKMASK  IS  TrUE  IN  THE  AREA  OF  *(,,*>  WHERE  COMPUTATION  TAKES  PLACE 
MASKMASK  ■  ROWS(2,YpnlNTS*l )  .AND,  COLS! 2, XpolNT$-t) 

T  ■  .FALSE, 

T  ■  ROWS! YRPOINTS,  YPOJNTS)  ,ANd.  COLS( XFPOJnTS,  xPoinTS) 
mASKMASK(T)  ■  .FALSE. 


MASKRB  IS  THE  REO/BIACK  SCHEME 

MASKRB  ■  ALTC(l)  .LEO,  ALTR! 1 ) 
MASKRB!. NOT.  MASKMASK)  «  .FALSE, 

TRACE  127  ( MERGE! 1 ,O.TRAN(MASKMASK )) ) 
TRACE  127  (MERGE! 1,0, TRAN(MASKRB))) 
RETURN 
END 


o  o  r» 


C  INITIALIZE  THE  »  MATRIX  - -  -  . 

c 

SUBROUTINE  INlT  W  :  _  •.  .r  -'g-  : 

COMMON  /ISCA/  ZPLANES,  XPOINTS,  YPOlNTS,  XFPOINTS,  XBPOINTS 
COMMON  /ISCA/  YRPOINTS, .YLPOINTS,  OAMHEISMT*.  0*mpace,  OAMRSIOE 
COMMON  /ISCA/  DAMLBACK,  DAMRBACK,  DAMLFSIOE,  OAMLBSIOE 
_  _  COMMON  /ISCA/  OAMMAXITERS,  BOTTOMMAXITERS,  NUMB I TERS,  NOMBOT  _ 

INTEGER  ZPLANESi  XPOINTS,  YPOINTSi  XFPOINTS,  XBPOINTS 
INTEGER  YRPOINTS,  YLPOINTS*  DAMHEIGHT,  0AMFAC6,  OAMRSIOE  __ 
INTEGER  OAMLBACK,  DAMRBACK,  DAMLFSIOE,  OAMLBSIOE  _ 

INTEGER  OAMMAXITERS,  BOTTQMMAXITERS,  NUMBITERS  ‘ 

PEAL  DAMEPSILON,  BOTTOMEPSILON,  OMEGA,  MAXDIFF 

COMMON  /RSCA/  DAMEPSILON,  BOTTOMEPSILON,  OMESA , lMAXOI FF  _ 

COMMON  /RMAT/M* , , 29 ) 

COMMON  /SUBLMAT/MASKRBI,  ),  MASKMASK*,)  ^ 

LOGICAL  MASKRB ,  MASKMASK 

REAL  TEMPS,  TEMPS!  _ 

TEMPS  a  OAMHEIGHT  *  DAMHEIGHT  *  0,5 
00  10  X  ■  1,  ZPLANES 

W( , ,K )  a  0,0 

TEMPS1  a  1  -  <K  -  1)  /  EFLOATI DAMHEIGHT )  _ 

W( ,1,K>  a  TEMPS  *  TEMPS1  *  TEMPS1 
10  CONTINUE 
RETURN 

ENO  _  ... 

INITIALIZE  the  bottom  z  PLANE  (I, E.  M(m1M 


SUBROUTINE  init  bottom  matrix 

COMMON  /ISCA/  ZPLANES,  XPOINTS,  YPOlNTS,  XFPOINTS,  XBPOINTS 
COMMON  /ISCA/  YRPOINTS,  YLPOINTS,  DAMHEIGHT,  DAMFACE,  OAMRSIOE 
COMMON  /ISCA/  DAMLBACK,  DAMRBACK,  DAMLFSIOE,  OAMLBSIOE 
COMMON  /ISCA/  OAMMAXITERS,  BOTTOMMAX I TERS ,  NUMBITERS,  NUMBOT 
INTEGER  ZPLANES,  XPOINTS,  YPOlNTS,  XFPOINTS,  XBPOINTS 
INTEGER  YRPOINTS,  YLPOINTS,  DAMHEIGHT,  DAMFACE,  OAMRSIOE 
INTEGER  OAMLBACK,  DAMRBACK,  DAMLFSIOE,  OAMLBSIOE 
INTEGER  OAMMAXITERS,  BOTTOMMAX ITERS,  NUMBITERS 
REAL  DAMEPSILON,  BOTTOMEPSILON,  OMEGA,  MAXOIFF 
COMMON  /RSCA/  DAMEPSILON,  BOTTOMEPSILON,  TMEGA,  MAXOIFF 
COMMON  /RMAT/W< , ,25) 

COMMON  /SUBlMAT/MASKPBI ,  MASKMASK*,) 

LOGICAL  MASKRB,  MASKMASK 
REAL  Z(.),  SAVEW* ,  )  *  ALPHA,  BETA 
LOGICAL  DONE 
INTEGER  NUMBTIMES 
NUMBOT  a  0 

ALPHA  a  OMEGA  *  0.25 
BETA  a  1.0  -  OMEGA 
10  SAVEW  ■  W( , ,  1  ) 

NUMBOT  a  NUMB  OT  ♦  1 
00  45  NUMBTIMES  a  1,2 
MASKRB(MASKMASK)  a  .NOT,  MASKRB 
a  W ( 3 , ,  i  ) 

Ml YPOlNTS, ,  1  >  ■  H( YP0INTS-2, , 1  ) 

Z  ■  W( , ,  1)  ♦  Ml-,-, 1  ) 

Z  a  Z(  +  ,  )  ♦  Z< »♦) 

Ml  MASkRr, 1 )  a  ALPHA  *  Z  ♦  BETA  *  H(,,l) 


i 
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A9  CONTINUE  _ 

maxdjff  •  MAiiAftiitAvew  -  w(##n)> 

DONE* ( MAXD IFF ,LE,B0TT0mEPS1 LON ) , OR , ( NUMBOT .GT.&UTTUmmAXITERS) 
If  (.NOT.  DONfi  «0TO  40 
TRACE  1 27 ( MAXOIFF  *  NUMBOT) 

return  _ _ _ 

END 

The  MAIN  LOOP  -  PROCESS  ALLLTHI  *  PLANES 
SUBROUTINE  HA  IN  LOOP 

COMMON  /ISCA/  ZPLANESi  XPOInTS#  YPOlNTS#  XPPSUnTS*  XBP01NTS 
COMMON  /ISCA/  YRPOINTS,  YLPOINTS,  OAMHglGHT,  DAMFACE#  oamrside 
COMMON  /ISCA/.  OaMLBAC^ a  OAMRftACK,  OaMLFSIDE#  OAML§SiOe  .. 

__  COMMON  /ISCA/  DAMMAXITERSl.  &OTTOMMAX ITERSi  NUMBlTERSf  NUMBOT 
INTEGER  ZPL ANES*  XPOINTS#  YPTSM2,  YPOlNTS#  XpPOlNTS#  XBPOINTS 
INTEGER  YRPOINTS#  YLPOINTS#  PAMHEICHT#  DAMFJC6#  OAMRSIDE 
1NTE6CR  OAMLBACK,  DAMRBACK,  DANLFSIDE#  DAMLBSIDt 
INTEGER  DAMMAXITERS,  BOTTOMMA X  ITERS ,  NUMBITERS 
real  OAMEPSILON,  BOTTOMEPSJLON#  OMEGA,  maxoipp 
COMMON  /RSCA/  OAMEPSILON#  BOTTOMEPSILON#  OMEQA#  MAXOIFF 

COMMON  /RMAT/NI  #»2S)  .  .  ."  _ 

COMMON  /SUBlMAT/MASKRBI , )#  MASHAS*!,) 

LOGICAL  MASKRBi  MASXMASK 
REAL  SAVEW( #  )#  Z( »  )#  Z|(« ) 

REAL  MAXSOFAR#  ALPHA#  beta#  W10GRI02#  WIDTHGrID 
INTEGER  NUMBTIMES#  TQPPLANE 
_  LOGICAL  TEMpMASKIt)#  DONE#  MSlQN(f) 

EQUIVALENCE  (NSiGNf Z) 

ALPHA  ■  OMEGA  *  1.0  /  6.0 
BETA  *  1.0  *  OMEGA 

WIDTH  OF  GRID  (I,E.  ONE  UNIT  SQUARE)  IS  SET  TO  1,0 

NUMBITERS  ■  0 

WIDTHGRID  ■  1,0 

WIDGRID2  •  WIOTHGRIO  *  WIOTHGRlD 
VPTSM2  ■  YPOlNTS  «.  2 
TOPPLANg  •  ZPLANES  -  1 

SAVE  THE  MASKRB  FOR  LATER  RESTORATION 

* 

TEMPMASK  ■  maskrb 

MaXQI.'F  is  the  MAXIMUM  DIFFERENCE  BETWEEN  SAVEW  AND  W(##K)  AFTER 
W(#.K>  HAS  ITS  RED  (OR  BLACK)  VALUES  CHANGE  (FOR  ALL  K) 

MAXOIFF  *  0.0 
NUMQITERS  ■  NUMBITERS  ♦  1 
MASKRB  a  TEMPMASK 
00  50  NJM8TIMES  ■  1,2 

ITERATE  FROM  THE  2N0  PLANE  TO  THE  TOP  PLANE 

DO  20  K  *  2#  TOPPLANF 
SAVEM  ■  «t(  »  ,K  ) 

REVERSE  RED/BLACK  FOR  SUCCESSIVE  PLANES 
MASKRB(MASKMASK)  ■  .NOT,  MASKRB 
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c 

e  SUM  TMf  III  MC!«M«0*S 
C 

$AV£»(lf)  ■  SHN(SAVEMlf) 

<?AVEM(  YPOINTS#  >  ■  SHS<SAVEMi2) 

1  ■  SAVEttt»#«!  s.^..  _ 

Zl  ■  SkyEM 

Zl(MASKRB)  ■  NCi, KM  ) 

Z(MASKRB)  ■  M( , ,KM ) 

^  |  •  *  ±  li  .Wv-Tt--  _ *_  ^ 

Zi  •  lit*)  .  _  _ 

. ,  _ «f. NOT. MASKRB)  *  .HlOfftlO*  x .  i  _ 

z  ■  z  ♦  zt 

Ij  !  *  Z(*» ) 

c 

ft  _  STORE  THE  ftVERAGE  OF  THE  SIX  NEIGHBOR*. IN  M  QNjL»  IN  THE 
ft  REO  (OR  BLACK)  CELLS 

z  Vhpha  *'  2  ♦  BETV  •"savgw  “ 

Z(NSISN)  ■  0.0 
M ( MASKRB  * K )  ■  Z 

c  .  . 

C  FIND  THE  MAXIMUM  DIFFERENCE  ON  THIS  PLANE 

e 

MAXSOFAR  ■  MAX(ABS(SAVEH  •  Z), MASKRB) 

IF  (MAXSOFAR  ,QT,  MAXOIFF)  MAXOIFF  •  MAXSOFAR 
20  CONTINUE 

ft.  _ _ 

C  REVERSE  STATE  OF  ORIGINAL  MASKRB  FOR  THE  2ND  PASS  THROUGH 

ft  the  planes 

C 

MASKRB(MASKMASK)  ■  .not,  tempmask 
30  CONTINUE 

DONE  •  (NUM3ITERS, AT. DAMMAXITERS), OR, (MAXOIFF, LE.DAMEPSILON) 

IF  (.NOT.  DONE)  GOTO  10 

MASKRB  ■  TEMPMASK 

RETURN 

END 
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APPENDIX  D:  The  three-dimensional  dam  problem  -  implementation  2. 


MASTER  EC ALPAL 

e 

C 

COMMON  /RMAT/MI 31, Jf) 

COMMON  /1SCA/  ZPLANESi  XPOINTS,  YPOINTS,  XFPOlNTf,  I»P0INT| 
COMMON  /!SCA/  TRPOINTS,  TLPOINTS,  OAMMfJCHf,  OA*FACEi  OAMtlOC 
COMMON  /ISC A/  0AML8AC*,  DAMRBACK,  DamLPSIDC,  DAMLBSIOE 
COMMON  /ISCA/  OAMMAXITCRSi  ftOTTOMMAI lT|RS«  NUtf»ITfRS#  NUMfOT 
INTEGER  2PLANES,  XPOINTS,  YPOINTS,  XPPOINTS,  XRPOINTS 
INTEGER  VRPOINTS,  TLPOINTS,  DAMME  I  GMT*  OAMfftCf*  OAMRSlQi 
INTEGER  OAMLGACK,  DaM»6AC*,  DAMLPSIOEi  DAMLBSIOE 
INTEGER  OAMMaXITERSi  BOTTOMMAX I T|R* ,  NUMRITERSi  TIMES 
PEAL  DAMEPSILON,  BOTTOMEPS ILON,  OMEGA,  MAXOIPP 
COMMON  /RSCA/  DAMIPS1LON,  BOTTOMfPSlLON,  OMCBAi  MAlDJFP 
INTEGER  V I ( 32  ) 

REAL  VRI 32 ) 

EQUIVALENCE  ( V  I ( ! ) , ZPL ANES ) , ( VR ( 1 ) , OAMEPS1LON) 

PAUSE  99 

10  REA0(?,100)  DAMHEIGMT,  DAMLPSIOEi  DAMLBSIOE.  DAMLBACK,  DAMRBaCK 

IPIDAMMEIGHT  .LE.  0)  PAUSf  00 
REA0(2,110)  MAXITERS.  EPSILON#  TIMES 
OAMEPSILON  •  EPSILON 
BOTTOMEPSILON  ■  EPSILON 
OAMMAIITERS  ■  MAXITERS 
BOTTOMMAXITERS  ■  MAXITERS 
NUMBIT1RS  •  1 
OMEGA  ■  1.6 

ZPLANES  ■  OAMHEIGHT  ♦  1 
OAMFACE  ■  DAMLBACK  ♦  DAMRBACK 
OAMRSIDE  ■  DAMLFS10E  ♦  DAMLBSIOE 
XPOINTS  «  OAMRSIDE  ♦  1 
XFPOINTS  *  DAMLFS IDE  ♦  1 
XBPOINTS  ■  DAMLBSIOE 
TRPOINTS  ■  DAMRBACK  ♦  2 
TLPOINTS  ■  OAMLBACKal 
XPOINTS  ■  TRPOINTS  ♦  TLPOINTS 

MR  I TE( 6# 1 20 )  DAMHEIGMT,  DAMLFSlDE#  DAMLBSIOE,  DAMLBaC*,  DamrbaC< 
CALL  OAPTO  IV( VI ,4hISCA,0) 

CALL  DAPTO  ES<DAMEPSILON,4MRSCA,0» 

CALL  DAPTO  ESI BOTTOMEPSILON, 4HRSCA# 1) 

CALL  DAPTO  ESI  OMEGA , 4HRSCA , 2 )  - 

CALL  OAPTO  IS<TIMES,4HR$CA,4) 

THF  TIME  OUT  PERIOD  IS  SET  TO  600  SECONDS 

CALL  DAPGOI 7MLAPLACE,  6001 
CALL  DApFROm  EV(VR,4MR$CA,0) 

CALL  OAPFROM  I VI V I i 4HISCA , 0 ) 
wR ITFI6, I  30)  NUMBITERS,  NUMBOT,  MAXDIPF 
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c 

C  rfPtTE  OUT  THE  M  MATRI*.  FROM  TOP  TO  BOTTOM  PUNE 
C 

NA  «  ZPLANES 
MI  •  VPOINTS 
NJ  •  APOIMTS 
DO  40  K  a  1 , NX 
-»ITE  (6.150)  4 

CALL  OAPPROM  EM(M,4HRMAT.32*(K«n  ) 

00  40  I  ■  1 , n 1 

WRITE  (6.140)  ( M( I i J ) i Jal « NJ ) 

40  CONTINUE 
GOTO  10 

00  50  J  a  1,NJ 
-RITE  (6.16c)  J 
DO  50  K*  ■  1  ,  NA 
<  •  NA«AA«1 

CALL  DAPFROm  EM(  W, 4HRMAT , 32*( A«1  )  ) 

WRITE  (6.140)  <W( I(J)iIaliNI) 

50  CONTINUE 

00  60  1  a  1.NI 
WRITE  ( 6 . 1  TO )  1 
00  60  kkmI.NK 
*  aNK-KKM 

CALL  OAPPROM  EM( W. AHRMaT . 52*( r*i  )  ) 

WRITE  (6.140)  (W(I.J)i Jal.NJ) 

60  CONTINUE 

GOTO  10 

100  FORMAT! 510  I 

ttO  FORMAT! 10.  EO.O,  10) 

120  FORMAT (  9M1MEIGHT  a,I5| 

1 22m  QAM  LEFT  FRONT  SlOf  a,I5,20H0AM  LEFT  BACA  SIDE 
212M  LE*T  »4CA  a, IJ, 12HRIGHT  BACA  a, 15,//) 

MO  FORMAT! 23H  NUMBER  OF  ITERATIONS  a, 15. 

19H  NUMBOT  a, I5.20HMAMIMUM  DIFFERENCE  a,E13,6) 

110  FOBMATdM  «  9E 1  3 . 6  ) 

160  FORMAT! ///I  IN  FOR  PLANE  ,|3> 

170  FORMAT! ///I 1M  FOR  SIDE  ,13) 

140  FORM AT ( ///l 1H  FOR  FACS  ,13) 

F  NO 
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swuesnMt.  lrpuce 

COMMON  /ISCA/  ZPLANES,  XPOInTS,  YPOInTS,  XFPOINTS,  XBPOINTS 
common  /ISCM  YRPOINTSr  YLPOINTS,  OAMHEUHT,  OAMFACE.  DAMRSIDE 
common  /ISC A/  DAMLBACK  ,  OAMRBACK,  damlfsioe,  oamlbsioe 
COMMON  /I$CA/  OAMMAXITERS.  BOTTOMMAX I TEPSi  NUMSlTERSr  NUMBOT 
INTEGER  ZPLANES,  XPOtNTS,  YPOiNTSi  xfpqints,  xbpoints 
INTEOifr VRRgINTS,  YLPOINTS,  OAMHEIGHT,  OAMFACE,  DAMRSIDE 
INTEGER  DAMLBACK,  DAMRBaCK,  DAMLFSIOE,  OAMLBSIOE 
INTEGER  OAMMAXITERS.  BOTTOMMA* ITE»$,  NUMftlTERS.  TIMES 
REAL  OAMEPSILON,  BOTTOMEPS tLON .  OMEGA.  MAXOIFF,  Z( ,  ) 

COMMON  /RSCA/  OAMEPSILON,  80Tf0MEPS_ILON,  OMEGA,  MAXOIFF,  TIMES 
COMMON  /RMAT/MI , ,25) 

COMMON  /SUdLttAT/MASXRSI , I,  MASKMASK! , ) 

LOGICAL  MASKRB,  MASK MASK  _ 

CALL  INIT  MASK 

CALL  INIT  N  _  _  _  .  , _ =  _ 

CALL  INIT  BOTTOM  PLANE 
CALL  MAIN  LOOP 

DO  1001  K«_2.ZPLANES.2,2_  _ 

2  ■  Wf»,Kf  _  .  _ _ _ 

WIMASKRB.K) »_W<  ,,kmi  __ 
w ( MASKRB, *>i >  •  Z 
1001  CONTINUE 
RETURN 
END 

C  INITIALIZE  the  masks 

SUBROUTINE  INIT  MASK 

COMMON  /ISCA/  ZPLANES,  XPOINTS,  YPOInTS,  xfpoints,  xbpoints 
COMMON  /ISCA/  yrpoints,  ylpoints,  damheigmt,  oamface,  damrside 

common  /ISCA/  DAMLBACK,  DAMRBACK,  DamLFSJDE,  OAMLBSIOE 
COMMON  /ISCA/  OAMMAXITERS,  BOTTOMMAX |TERS,  NUMBITERS,  NUMBOT 
INTEGER  ZPLANES,  XPOINTS,  YPOlNTS,  XFPOINTS,  XBPOINTS 
INTEGER  YRPOINTS,  YLPOINTS,  DAMHEIGMT,  OAMFaCE,  DAMRSIDE 
INTEGER  DAMLBACK,  DAMRBACK,  DAMLFSIOE,  OAMLBSIOE 
INTEGER  OAMMAXITERS,  BOTTOMMAX ITERS,  NUMglTEPS,  TIMES 
REAL  OAMEPSILON,  BOTTOMEPSILON,  OMEGA,  MAXOIFF 
COMMON  /RSCA/  OAMEPSILON,  BOTTOMEPSILON,  OMEGA,  MAXOIFF,  TIMES 
COMMON  /RMAT/Wt , ,25) 

COMMON  /SUlLMAT/MASKPRJ , ) I  MASKMASK ( , ) 

LOGICAL  MASKRB,  MASKMASK 
LOGICAL  T( , ) 

MASKMASK  IS  TRUE  IN  TME  area  OF  W(,,K)  WHERE  COMPUTATION  TAKES  PLACE 


v.  a  S  <  MASK  ■  ROWS!  2,  YPOINTS-1)  .AND,  COLS<  2,  XPPINTS-1  ) 

T  ■  RO*S< YRPOINTS,  YPOInTS)  .AND.  COL S ( XPPO l NTS ,  XPOINTS) 
MASkMASK(T)  »  .FALSE, 


MASKRR  IS  THE  REO/RLAC*  SCMEMp 

MASKRB  «  ALTCU)  .LEO.  ALTRM) 

IF  ( . NOT. SWITCH! 2 ) )  MASpRB* . NOT . MASKRB 
TRACE  127  ( MFRGEI 1 , :,TWANIMAS*MASK  )  )  J 
TRACE  127  IMERGE(1,?,TRAN(MASKRR)») 

return 

END 
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u  IMTULIZC  Tw6  «  MA.-UX 
C  .  ...... 

.  $UBRQO?IN«  I  NIT  * 

COMMON  /ISCA/  ZPLANES,  XP0INTS,  YPOInTS,  APPOINTS,  XBPOINTS 
COMMON  /ISCA/  APPOINTS#  YLPOINTS,  DAMHEIGHT,  OAMFACE,  OAMRSIDE 
COMMON  /ISCA/  DAMLBACK ,  DAMRBACK ,  OamlFSIOE,  DAMLBSIOE 
COMMON  /ISCA/  DAMMAXITERS,  BOTTOMMAXtTERS,  NUMBITERS, __NUMBOT 
- ‘  EGER  ZPLANES,  XPOlNTS,  YPOINTS,  XF>OTnTS,  XBPOINTS 
IN  EGER  YRPOINTS,  YLPOINTS,  OAKHEISHT*  DAMFACEi  OAMRSIDE 
DAMLBACK i  DAMRBACK ,  OAMLFSIOE,  OAML8S I OE 

INTEGER  dammaxjters,  bottommaxiters,  numbiters,  times 

real  OAMEPSILON,  BOTTOMEPSILON,  OMEGA ,  MAXDIFF,  Z( ,  ) 

COMMON  ;  /R  $c  A/  ftAMf PS  I L  ON ,  BOTTOMEPSILON,  OMEGA,  MAXDIFF,  TIMES 
COMMON  /RMAT/MI , ,29) 

COMMON  /SOBLMAT/MASkRBI  , ),  MASKMASK!,) 

LOGICAL  MASKflB,  MASKMASK 

rial  temps,  tempsi 

TEMPS  ■  DAMHElSHt,*  DAMHEIGHT  *  0,9 
DO  JO  X  ■  1,  ZPLANES 
2L  *<#.K)  *  0.0 

TEMPSI  ■  1  -  (K  ,  n  /  EFLOaTIDAMHEIGHT) 

Ml*l,K)  ■  TEMPS  *  TEMPSI  *  TEMPSI 
10  continue 

00  tool  Ks  2, ZPLANES»2, 2 

Z  ■  Ml , ,K ) 

i  WtMASMRB,K)  m  W|,,K*1) 

M(MASKRB,K*1 )  ■  Z 

iOoi  continue 

W(  •  .ZPLANESM  )  ■  M(,,  ZPLANES) 

RETURN 

END 


INITIALIZE  THE  BOTTOM  Z  PLANE  II. E.  M|,,D) 

.subroutine  INIT  bottom  matrix 

COMMON  /ISCA/  ZPLANES,  XPOInTS,  YPOInTS,  appoints,  xbpoints 
COMMON  /ISCA/  YRPOINTS,  YLPOINTS,  DAMHEIGHT,  OAMFACE,  OAMRSIDE 
COMMON  /ISCA/  DAMLBACK,  DAMRBACK,  OAMLFSIOE,  DAMLBSIOE 
COMMON  /ISCA/  DAMMAXITERS,  BOTTOMMA X I TERS ,  NUMBITERS,  NUMBOT 
INTEGER  ZPLANES,  XPOlNTS,  YPOINTS,  XfPQJNTS,  XBPOINTS 
INTEGER  YRPOINTS,  YLPOINTS,  DAMHEIGHT,  OAMFACE,  DAMHSIOE 
INTEGER  DAMLBACK,  DAMRBACK,  OAMLFSIOE,  DAMLBSIOE 
INTEGER  DAMMAXITERS,  BOTTOMmA X  I TERS ,  NUMBITERS,  TIMES 
REAL  OAMEPSILON,  BOTTOMEPSILON,  OMEGA,  MAXDIFF 

COMMON  /RSCA/  OAMEPSILON,  BOTTOMEPSILON,  OMEGA,  MAXDIFF,  TIMES 
COMMON  /RMAT/NI , ,29) 

COMMON  /SUBLMAT/MASkRBI , ) ,  maSKMASKI.) 

LOGICAL  MASKR8,  MASKMASk 

REAL  ?(.),  SAVEWI,),  ALPHA,  BtTA 

LOGICAL  OONE 

integer  numbtimes 
njmbot  ■  0 

ALPHA  ■  OMEGA  *  0.25 
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non  r>  nn  I 


W  :  .  . - - - - - - 1 


BET  ft  *  1.0  -  OMEGA 
SAVE*  ■  W  (  , ,  1  ) 
numbot  ■  NUMB  OT  +  i 
00  15  NUM1TIMES  »  1,2 
W(  i,,l)  ■  M( 5 , , 1 ) 

W(  YPOINTS,  •  1 )  ■  W(YP0INTS»2,,n 

2  ■  W< t , 1  )  ♦  W(*,w, 1  ) 

2  ■  2  C  ♦ »  )  ♦  2<  ,♦) 

W(MAS<RB. ANO.MASKMASX, 1)  ■  ALPHA  *  Z  ♦  BETA  *  W(,,l) 

MASKRB  ■  .NOT,  MASKRB 
CONTINUE 

MAXDIFF  •  MAXI ABStSAVEH  *  W(,#l))) 

DONE«<MAXDIFF,LE.BOTTOMEPSlLONJ .OP. (NUMBOT.GE, BOTTOMMAX ITERS) 
IF  ( .NOT.  DONE)  GOTO  40 
TRACE  1 27 ( MAX  0  IFF#  NUMBOT) 

RETURN 
END 

THE  MAIN  LOOP  -  PROCESS  All  THE  Z  planes 
SUBROUTINE  MAIN  LOOP 

COMMON  /ISCA/  ZPLAN6S,  XPOINTS,  YPOINTS,  XFP01NTS#  XBPOINTS 
COMMON  /ISC ft/  YRPOlNTSi  YLPOINTS.  DAMHEIGHT,  DAMFACE,  DAMRSIOE 
COMMON  /ISCA/  DAMLBACK #  OAMRBACK,  OAMLFSIOE.  DAMLBSIUE 
COMMON  /ISCA/  0AMMAX1TERS,  BOTTOMMAX I TERS ,  NUMBITERS,  NUMBOT 
INTEGER  ZPLANES,  XPOINTS,  YPTSM2,  YPOINTS,  XFPOJNTS,  XBPOINTS 
integer  YRPOINTS,  YLPOINTSi  OAMHEIGHTt  OAMFaCE,  DAMRSIOE 
INTEGER  DAMLBACK,  OAMRBaCK,  OAMLFSIDB#  DAMLB5I0E 
INTEGER  OAMMAX I TERS,  BOTTOMMAX ITERS,  NUMBITERS,  TIMES 

peal  damepsilon,  bottomepsilon,  omega#  maxdiff 

COMMON  /RSCA/  DAMEPSILON,  BOTTOMEPSILON,  OMEGA,  MAXDIFF,  TIMES 
COMMON  /RMAT/W( . ,25) 

COMMON  /SUBLMAT/MASKRBI ,  ) #  MASK  MASK ( ,  ) 

COMMON  /WORK/  Z,WK 
LOGICAL  MASKRB,  MASKMASK 

REAL  SAVEW(,),  Z(,),  ZKili  WKPJ(  ,  ),  WK  (  i  )  ,  MAXO(  ,  ) 

REAL  MAXSOFAR,  ALPHA,  BETA,  WIQGRID2,  WIqTHGRIO 
integer  NUMBTIMES,  topplane 

LOGICAL  TEMPMASKI,),  DONE,  WSIGNI,),  TEST,  NOTTEST 
EQUIVALENCE  (WSIGN.Z),  (ZiZl)#  (WKiMKPl) 

ALPHA  «  OMEGA  *  1,0  /  5, 0 
BETA  *1.0-  OMEGA 

WIDTH  nF  GRID  CI.E.  ONE  UNIT  SQUARE)  IS  SET  TO  1,0 

numriters  ■  0 

wIOfHGRIO  ■  1 ,0 

w I 0G° 102  ■  WIOTHGRIO  *  WIDTwGRID 
TOPPLANE  ■  ZPLANES  -  1 


40 


45 
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TEST  ■  TIMES  «EQ,  1 
HOTTEST  ■  .NOT,  TEST 
ITHES  ■  TIMES 
MAXD  ■  0,0 


ALTER  ALL  THE  000  NUBEREO  PLANESI 
K M2  «  1 

00  20  K  ■  2,T0PPLAN£,2 
SAVEW  ■ 

WK  ■  W( , #K  ) 

W<< i, J  ■  SHN( w«,2) 

WK  ( YPO I NTS # )  ■  SH$<WK.2) 

2  ■  wK  ♦  WK(»r*) 

7  ■  (?<♦, )  +  ll »*)*WK*MER6E<W( ,  .KM2)  iW( » »K  +  2) .maSKRB )-WlUGRI02)*ALPH 
Z(WSIGN)  ■  0.0 
W(MASKMASK*K+1 )  ■  2 
IP  (NOTTEST)  GOTO  20 
21  •  ABS(SAVEW«Z) 

MAXD(Zl.GT.MAXD)  ■  Z1 
0  KM2  ■  K 


ALTER  ALL  THE  EVEN  NUMBERED  PLANESI 
00  21  K  ■  2 , TOPPLANE  $  2 
SAVEW  ■  W(,,K) 

WKPl  ■  W<|'K+1) 

WKP 1 ( 1 f )  -  SHN(WKP1,2I 
WKPl ( VPOINTS# 1  ■  SH$( W*P1 • 2 ) 

2  ■  WKPl  ♦  WKPl(-i^) 

7  ■  {  7  (  ♦  *  )  *2  (  »*)  +  WKPl  ♦MERGE!  W(  »  ,03)  ,W<  #  ,K-1  ) .  MASK  RB  I w  I OGR  1 02  )  *  AL 
Z(WSISN)  m  0.0 
W(MASKMASKiK)  ■  Z 
IF  (NOTTEST)  SOTO  21 
71  «  ABS(SAVEW.Z) 

MAXD(71 .GT.MAXD)  ■  Z1 
1  CONTINUE 


ITHES  ■  ITIMES^l 
IF  ( ITIMES.GT, 1 )  GOTO  2 
IF  ( ITIMES.EO.O)  GOTO  3 
TEST  ■  .TRUE. 

NOTTEST  ■  .FALSE. 

GOTO  2 
C 

3  NJMBITERS  ■  NUMBITERS  ♦  TIMES 

MAXOIFF  *  MA¥(MAXO»MAS*MASK> 

IF  (MAXDIFF.LT.OAMEPSILON)  GOTO  4 
C 

IF  (NUMBITERS. LT.0AMMAX1TERS)  GOTO  l 

4  RETURN 

FNO 
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